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ABSTRACT 


The localization of mobile wireless communication units is studied. The most 
important method of localization is the time difference of arrival (TDOA) method. The 
wavelet transform is used to increase the accuracy of TDOA estimation. Several 
denoising techniques based on the wavelet transform are presented in this thesis. These 
techniques are applied to different types of test signals as well as the GSM signal. The 
results of the denoising techniques are compared to the ones obtained using no denoising. 
The denoising techniques allow a 28 to 81 percent improvement in the TDOA estimation. 
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1. INTRODUCTION 


The problem of providing reliable and accurate position information of mobile 
wireless communication units, has attracted a lot of attention in recent years. The main 
factor behind the recent interest has been the adoption of certain regulations by the 
Federal Communications Commission (FCC) [1], that require wireless communications 
licensees to incorporate a position localization capability in their systems to provide 
Enhanced-911 (E-911) service. However, there are many other reasons for wireless 
service providers to have such a system in place. For example, one can use reliable 
position localization as means to optimize the performance and design of the wireless 
networks. 

A. THESIS OUTLINE 

In the remainder of this Chapter, reasons for localizing cellular emitters and 
standard methods to accomplish this task are explored. Chapter 11 introduces the GSM 
system. Chapters III-VI discuss wavelet based denoising methods. Test signal description 
and simulation results are contained in Chapter Vn. Chapter VIII contains the conclusion 
and suggests logical extension of the work. 

B. THESIS CONTRIBUTIONS 

The main contribution of this work is the reduction in mean square error for the 
time difference of arrival estimate. This permits localization with corresponding smaller 
error. The novel denoising methods, a modified approximate maximum likelihood 
(MAML), a fourth order moment, and a time varying MAML method, provide desirable 
improvements. 
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C. WHY IS THERE A NEED TO LOCATE CELLULAR PHONES? 

1. Public Safety and Enhanced Emergency Services 

The FCC of the United States passed a regulation to provide E-911 service 
whereby “carriers are required to identify the latitude and longitude of a mobile unit 
making a 911 call.” At present, if a cellular subscriber dials 911 and does not relay 
location information to the operator, the authorities can only estimate the caller’s location 
to within a few kilometers, based on the cell being used. Providing location and tracking 
to emergency services would save critical seconds in response time to stranded motorists, 
injury victims/witnesses who are confused or unable to relay geographical information, or 
allow vehicle pursuit by authorities. 

2. Fleet/Asset Management for Couriers and Transportation Businesses 

Many organizations already utilize cellular phones for their day-to-day business 
activities. It makes sense to utilize these phones, or wireless transmitters using cellular 
technology, to monitor and track vehicles and shipments, including couriers, taxis, 
buses, and other fleet-based commercial services. 

3. Tracking of Stolen Phones and/or Criminals 

Stolen and cloned cellular and personal communications systems (PCS) phones 
represent a major problem for cellular carriers and users. Authorities acknowledge an 
increased use of mobile phones in criminal activities. Until now, mobile phones used in 
illegal activities have been difficult to trace and millions of dollars in wireless fraud 
have been committed. With the availability of cellular phone tracking and localization, 
the criminal element in wireless phone use could diminish dramatically. 
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4. Tracking Stolen Vehicles 

Many vehicles are equipped with a cellular phone, and many vehicles are being 
manufactured with a tracking “tag” for location or navigational purposes. The 
localization technology and its future variations offer the ability of recovering lost or 
stolen vehicles. 

5. Location-Sensitive Navigation 

Imagine driving into new and unfamiliar city, and being able to dial up 
navigational service. The operator could identify your current location and give you 
directions to the nearest service station or hotel based on your present position. 

6. Localization of Cellular and PCS Telephone Users 

Similarly to the navigational service, a pay-per-use service may allow tracking of 
persons who use cell phones. Parents who need to locate their children or family 
members wanting to find a medical patient will benefit from this technology, saving time 
and reduce worry. If one is uncertain about his own location, this service could help him 
to give position or address information. 

7. Electronic Databases 

Often, it is useful to study the demographics of an area to determine the necessity 
for roads and other infrastructures. By tracking the use and location of cellular phone 
users, cellular carriers can strategically plan base stations. 

8. Law Enforcement or Military Use 

Position information can also be made available to law enforcement and military 
users. The information allows localization of a user independent of his desire to do so. 
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D. EXISTING POSITION LOCALIZATION SYSTEMS 


A number of position localization systems have evolved over the years. These 
position localization systems are: Global Positioning System (GPS), Loran C, SIGNPOST 
Navigation, Global Navigation Satellite System, Automatic Vehicle Monitoring, and 
Cellular Phones Geolocation. 

1. Global Positioning System (GPS) 

GPS is the most popular radio navigation aide, due to high accuracy, worldwide 
availability, and low cost. GPS is also used to relay position of cellular phones to public 
switched telephone networks (PSTN) or to public safety answering points (PSAP). 

The principle behind GPS is very simple. GPS uses a time-of-amval (TOA) 
method. GPS uses precise timing within a group of satellites that transmit a spread 
spectrum L-band signal to ground in the centered at 1575.42 MHz. An accurate clock at 
the receiver measures the time delay between the signals leaving the satellite and arriving 
the receiver. This allows calculation of the distance from the satellite to the cellular 
phone. If one obtains three distances using three different satellites, one can use a 
triangulation method to determine the position of the cellular phone. 

Currently a GPS receiver costs less than $200, and its accurate to 100 meters [13, 
23]. More sophisticated units, including those used by military, which use differential 
GPS, provide an accuracy of a few meters. 

2. Loran C 

Loran C is a navigational tool that operates in the low frequency (90-110 kHz) 
band and uses a pulsed hyperbolic system for triangulation. It has a repeatable accuracy 
in the 19-90 meters range and is accurate to about 100 meters with 95 percent confidence 
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and 97 percent availability. like GPS, its performance depends on local calibration and 
topography. GPS has replaced Loran C in most applications [13]. 

3. SIGNPOST Navigation 

SIGNPOST Navigation employs a large number of simple radio transmitters to 
accurately determine position of a mobile. These transmitters are replaced along 
highways and typically serve as coded beacons, where the code designates the latitude 
and longitude of the SIGNPOST. The transmitter strength indicates the relative position 
of the receiver to the transmitter. This navigation aids work well for limited areas such as 
a small city. While not originally designated as such, today’s Advanced Mobile Phone 
System (AMPS) analog cellular radio system may actually serve as a SIGNPOST system, 
since each base system transmits a beacon signal on its forward control channel [24]. As 
a part of the forward control channel structure, an overhead message containing a station 
identification number (SID) and a digital color code (DCC) is sent every 0.8 seconds. 
The DCC is used to determine the location within the cell [13]. 

4. Global Navigation Satellite System (GLONASS) 

The Global Navigation Satellite System (GLONASS) is similar system to GPS. 
Although the system uses principles similar to GPS, its operation differs in several 
aspects. The synchronization period of GLONASS takes only 1/3 as long as GPS, 
typically under a minute [13]. 

5. Automatic Vehicle Monitoring (AVM) 

Automatic Vehicle Monitoring (AVM) systems provide position localization 
capabilities for handling a large number of vehicles. Typical applications include fleet 
management, vehicle security, and emergency services. AVM systems have been 
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available in the United States for a number of years [13]. In 1995, the FCC changed the 
name of these systems to Location and Monitoring Services (LMS). In the United States, 
the primary band for LMS is the 902-928 MHz industrial, scientific, and medical (ISM) 
band. LSM is also supported to a lesser extent in several bands below 512 MHz. The 
LSM system is a licensed system with up to 300 W peak power for the forward link; 
however it shares the band with low-power unlicensed devices, such as cordless phones, 
wireless local area networks, and utility meter reading systems. The band is also used by 
federal government radiolocation systems and amateur radio operators, so the prospect of 
the interference between LMS and other users of the spectrum is an issue in the 
deployment of LMS systems [25]. 

E. METHODS FOR LOCATING CELLULAR PHONES 

There are many techniques that can be considered for the localization of cellular 
phones. These techniques can be classified into two categories. These are, position 
localization systems that require a modification of the existing handsets, and systems that 
require modification at the base stations. The modification of existing handsets can be 
accomplished by using the GPS based position localization system, mobile assisted time 
or time difference of arrival techniques. The second category consists of angle of arrival 
(AOA), time difference of arrival (TDOA), time of arrival (TOA), or frequency 
difference of arrival (FDOA) estimation of the wavefront at the receiving platforms. 

1. Angle of Arrival (AOA) 

AOA can be accomplished by means of a highly directive receiving antenna or by 
means of a nulling measurement using several feeds from the antenna. A single platform 
may be sufficient for AOA position localization of a wireless transmitter on the surface of 
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the earth. Once the angle is precisely determined, the position of the emitter can be 
determined by the intersection of the centerline of the antenna beam, the boresight, with 
the surface of the earth. 

Although AOA methods offer a practical solution for the wireless transmitter 
localization, they have certain drawbacks. For accurate AOA estimates, it is crucial that 
signals coming from the source to antenna arrays must be coming from the Line-Of-Sight 
(LOS) direction. However, this is often not the case in cellular systems, which may be 
operating in heavily shadowed channels, such as those encountered in urban 
environments. Another factor is the considerable cost of installing antenna arrays. This 
method also requires a relatively complex AOA algorithm. Although there are 
exceptions, these algorithms tend to be highly complex because of the need for 
measurement, storage and usage of array calibration data and their computationally 
intensive nature [3,13]. 

2. Frequency Difference of Arrival (FDOA) 

FDOA measurements require at least two receiving platforms and that the relative 
velocity between the platforms be large enough that the difference in Doppler shifts of 
the two received signals is significantly larger than the frequency measurement error. 

3, Time of Arrival (TOA) 

It may be possible for the base station to indirectly determine the time that the 
signal takes from the source to the receiver on the forward or the reverse link. This may 
be done by measuring the time in which the mobile responds to an inquiry or an 
instruction transmitted to the mobile from the base station. The total time elapsed from 
the instant the command is transmitted to the instant the mobile response is detected, is 


7 



composed of the sum of the round trip signal delay and any processing and response 
delay within the mobile unit. If the processing delay of the response within the mobile is 
known with sufficient accuracy, it can be subtracted from the total measured time, which 
provides us the total round trip delay. 

There are certain problems with this method. The estimate of the response delay 
within the mobile might be difficult to determine in practice. The main reason is the 
variations in designs of the handsets from different manufacturers. Secondly, this method 
is highly susceptible to timing errors in the absence of LOS, as there would be no simple 
way to reduce the errors induced by multiple signal reflections on the forward or reverse 
link. 

4. Time Difference of Arrival (TDOA) 

The classical approach to estimating TDOA is to compute the cross correlation 
between signals arriving at two base stations. The TDOA estimate is taken as the delay, 
which maximizes the cross-correlation function. The cross-correlation function is also 
used to determine at which base station the signal arrives first. These two pieces of 
information yield a hyperbolic localization curve. We can localize the wireless 
transmitter by solving two hyperbolic curve equations. 

It is necessary that the code generators at each receiver be synchronized so that 
the TDOA estimates have a common time base [2]. This form of radio localization is 
useful in asynchronous system since time of transmission need not be known. In 
geometric interpretation, this procedure reduces to finding the intersection of hyperbolas 
whose foci are at the receivers. To determine the location of a transmitter in two 
dimensions, at least three receivers are required. 
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This method offers many advantages over other competing techniques. No 
modification of the existing handsets is required. In this respect, this solution would be 
more cost effective than a GPS-based solution. It also does not require knowledge of the 
absolute time of transmission from the handset like a modified TOA method needs. Since 
this technique does not require any special type of antennas, it is less expensive to put in 
place than the AOA methods. It can also provide some immunity against timing errors if 
the source of major signal reflections is near the mobile. If a major reflector effects the 
signal components going to all the receivers, the timing error may get cancelled or 
reduced in the time difference operation. Hence, TDOA methods may work accurately in 
some situations where there is no LOS signal component. In this respect, it is superior to 
the AOA method and TOA methods. 
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11. GLOBAL SYSTEM FOR MOBILE (GSM) 


GSM is a second-generation cellular system standard that was developed to solve 
the fragmentation problems of Europe’s first cellular systems. GSM is the first cellular 
system to specify digital modulation, network level architectures and services. Before 
GSM, European countries used different cellular standards throughout the continent, and 
it was not possible to use a given single subscriber unit throughout Europe. GSM was 
originally developed to serve as the Pan-European cellular service and promised a wide 
range of network service through the use of the Integrated Services Digital Network 
(ISDN). GSM’s success has exceeded the expectations of virtually everyone, and it is 
now the world’s most popular standard for new cellular radio and personal 
communications equipment. 

A. GSM SYSTEM ARCHITECTURE 

The GSM system architecture consists of three major interconnected subsystems 
that interact with themselves and the users through network interfaces. The subsystems 
are the Base Station Subsystem (BSS), Network and Switching Subsystem (NSS), and the 
Operation Support Subsystem (OSS). The Mobile Station (MS) is also a subsystem, but is 
usually considered to be part of the BSS for architectural purposes. 

The BSS provides and manages radio transmission between the MS’s and the 
Mobile Switching Center (MSC). The BSS also manages the radio interface between the 
MS’s and all other subsystems of GSM. 
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The NSS manages the switching functions of the system and allows the MSC’s to 
communicate with other networks such as the Public Switched Telephone Network 
(PSTN) and Integrated Services Digital Network (ISDN). 

The OSS supports the operation and the maintenance of GSM and allows system 
engineers to monitor, diagnose, and troubleshoot all aspects of GSM. This subsystem 
interacts with the other GSM subsystems. 

1. GSM Signal SpeciHcations 

GSM utilizes two 25 MHz bands, which have been set aside for system use in all 
member countries. The 890-915 MHz band is used for subscriber-to-base transmission 
(reverse link), and the 935-960 MHz band is used for base-to-subscriber transmission 
(forward link). GSM uses Frequency Division Duplex (FDD) and a combination of Time 
Division Multiple Access (TDMA) and Frequency Hopped Multiple Access (FHMA) 
schemes to provide stations with simultaneous access to multiple users. The available 
forward and reverse frequency bands are divided into 200 KHz channels. These channels 
are identified by their Absolute Radio Frequency Channel Number (ARFCN). The 
ARFCN denotes a forward and reverse channel pair, which is separated in frequency by 
45 MHz. Each channel is time shared between as many as eight subscribers using 
TDMA. 

The radio transmissions on both forward and reverse link are made at a channel 
data rate of 270.833 Kbps using binary 0.3 OMSK modulation. The bandwidth-bit 
duration product, BT, has a level of 0.3. The signaling bit duration is 3.692 p,s. User data 
is sent a maximum rate of 24.7 Kbps. Each time slot (TS) has an equivalent time 
allocation allowing for 156.25 channel bits. From this, 8.25 bits of guard time and 6 total 
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start and stop bits are used to prevent overlap with adjacent time slots. Each TS has a 
time duration of 576.92 ps, while a single GSM TDMA frame spans 4.615 ms. The total 
number of available channels within a 25 MHz bandwidth is 125. Table 2.1 summarizes 
the signal specifications. 


Parameter 

SpeciBcations 

Reverse Channel Frequency 

890-915 MHz 

Forward Channel Frequency 

935-960 MHz 

ARFCN Number 

0 to 124 and 975 to 1023 

Tx/Rx Frequency Spacing 

45 MHz 

Tx/Rx Time Slot Spacing 

3 Time slots 

Modulation Data Rate 

270.833333 Kbps 

Frame Period 

4.615 ms 

Users per Frame (Full Rate) 

8 

Time Slot Period 

576.9 ps 

Bit Period 

3.692 ps 

Modulation 

0.3 GMSK 

ARFCN Channel Spacing 

200 KHz 

Interleaving (max delay) 

40 ms 

Voice Coder Bit Rate 

13.4 Kbps. 


Table 2.1: GSM signal specifications. 


2. Gaussian Minimum Shift Keying (GMSK) 

GMSK is a binary modulation scheme which may be viewed as a derivative of 
Minimum Shift Keying (MSK). In GMSK, the sidelobe levels of the spectrum are reduced 
by passing the modulating Non-retum to zero (NRZ) data waveforms through a pre¬ 
modulation Gaussian pulse-shaping filter. Baseband Gaussian pulse shaping smoothes the 
phase trajectory of the MSK signal and hence stabilizes the instantaneous frequency 
variations over time. This has the effect of considerably reducing the sidelobe levels in 
the transmitted spectrum. 
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3. 


GSM Channel Types 


There are two types of GSM logical channels, called traffic channels (TCH) and 
control channels (CCH). Traffic channels carry digitally encoded user speech or user data 
and have identical functions and formats on both the forward and reverse link. Control 
channels carry signaling and synchronizing commands between the base station and the 
mobile station. 

4. Frame Structure for GSM 


As shown in Figure 2.1, there are eight time slots (TS) per GSM frame. The frame 
period is 4.615 ms. A TS consists of 148 bits which are transmitted at rate of 270.833 
Kbps (an unused guard time of 8.25 bit period is provided at the end of each burst). Out 
of the total 148 bits per TS, 114 are information bits, which are transmitted as two 57 bit 
sequences close to the beginning and end of the burst. A 26 bit training sequence allows 
the adaptive equalizer in the mobile or base station receiver to analyze the radio channel 



Figure 2.1: GSM frame structure. 
















characteristics before decoding the user data. Stealing flags on the both side of the 
training sequence are used to distinguish whether the TS contains voice or control data. 

B. TRANSMITTER/RECEIVER SYSTEM 

In this section we will look at the GSM transmitter and receiver and simulate 
some GSM signals to be used for simulation purposes. The overall structure of the GSM 
transmitter/receiver system is illustrated in Figure 2.2. 
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Figure 2.2:Block diagram for the GSM transmitter/receiver system. 


1. Transmitter Structure 

The overall stmcture of the transmitter is illustrated in Figure 2.3. The transmitter 
is made up for distinct functional blocks. 
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Figure 2.3: The overall structure of the transmitter. 


To simulate an input data stream to the channel encoder/interleaver, a sequence of 
random data bits is generated by the random bit generator. This sequence is accepted by 
the MUX, which splits the incoming sequence to form a GSM normal burst. The burst 
type requires that a training sequence is included which also must be supplied. Upon 
having generated the prescribed GSM normal burst data structure, the MUX sends this to 
the GMSK-modulator, where GMSK is short for Gaussian Minimum Shift Keying. The 
GMSK-modulator block perfonns a differential encoding of the incoming burst to form a 
NRZ sequence. This modified sequence is then subject to the actual GMSK-modulation 
after which, the resulting signal is represented as a complex baseband 1 and Q signal. 

2. Receiver Structure 

The general structure of the receiver, consisting of three functional blocks, is 
illustrated in Figure 2.4. The demodulator accepts the GSM burst, r, using a complex 
baseband representation. Based on this data sequence, the oversampling rate OSR, the 
training sequence TRAINING , and the desired length of the receiving filter, Lh, the 
demodulator determines a bit sequence. This demodulated sequence is then used as the 
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input to the demultiplexer (DeMUX) where the actual data bits are obtained. The 
remaining control bits and the training sequence are stripped off. Finally channel 
decoding and de-interleaving is performed. 



Figure 2.4: The overall structure of the receiver. 
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III. WAVELETS 


Usually signals are transformed to obtain information that is not directly 
observable in the raw signal. There are many transformations that could be used. The 
wavelet transform belongs to this set. It provides a time-scale representation [6]. There 
are other transformation, which can give similar information, such as the short time 
Fourier transform, Wigner-Ville distribution, etc. 

Wavelet analysis can be interpreted as an extension of the Fourier analysis. Thus, 
this chapter will first discuss Fourier and then present the wavelet analysis. 

A. FOURIER ANALYSIS 

Fourier analysis breaks a signal into its constituent sinusoidal components at 
different frequencies. One can also think of the Fourier analysis as a mathematical 
technique for transforming the signal from a time-based to a frequency-based 
representation [7]. 

1. Fourier Series 

Let gp(r)denote a periodic signal with period Tg. By using the Fourier series 

expansion of this signal, we can resolve this signal into an infinite sum of complex 
exponentials. The Fourier series expansion can be written as [8]: 

, (3.1) 

where n/Tg represents the nth harmonic of the fundamental frequency /g =HTq. The 
series expansion of Equation 3.1 is referred to as the complex exponential Fourier series. 
The c„ coefficients are the complex Fourier coefficients. We can determine c„ as follows: 
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(3.2) 


C. 



0 ,± 1 ,± 2 ,- 


The complex exponential analysis equation provides the coefficients necessary to 
reconstruct the periodic signal from its Fourier series coefficients. A plot of the 
magnitude of c„ versus frequency is called the magnitude spectrum of the signal g pit). 


The spectrum provides frequency information of the signal. 

2. Fourier Transform 

In the previous section we used the Fourier series to represent a periodic signal. 
We will develop a similar representation for a signal g(/)that is nonperiodic. In order to 
do this, we first construct a periodic function gpit) of period of T^, in such a way that g(t) 
defines one cycle of this periodic function. In limit we let the period Tq become infinitely 
large, so that we may write 

g(f)= limg (0 . (3.3) 

The Fourier transform of a general continuous function g(t) is defined as [8]; 


= . (3.4) 

G(f) is a continuous function of the frequency variable/. The original signal g(t) can be 
recovered exactly from G(f) by means of the inverse Fourier transform which is defined 
as follows [8]: 

= (3.5) 

The two functions g(t) and G(f) uniquely define each other and are known as a Fourier 
transform pair. 
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Fourier analysis is extremely useful because the signal’s frequency content is of 
great importance. So why do we need another techniques, such as the wavelet analysis? 

Fourier analysis has a serious drawback. In transforming to the frequency domain, 
time information is lost. When looking at a Fourier transform of a signal, it is impossible 
to tell when a particular event took place. If a signal is stationary, this drawback is not 
very important. However many interesting signals contain numerous non-stationary or 
transitory characteristics such as trends, abrupt changes, and beginnings and ends of 
events. These characteristics can be the most important part of the signal, and Fourier 
analysis is not suited to localize them. 

3. Short-Time Fourier Analysis 

The Fourier analysis technique described above provides a frequency domain 
presentation of the signal. When the signal is non-stationary, it is desirable to have a 
description that involves both time and frequency. 

The short-time Fourier transforms (STFT) can be viewed as an extension of the 
Fourier transform devised to map the signal into the two dimensional time-frequency 
plane. The STFT uses a sliding window function M>(t) to segment the signal into small 
uniform blocks of time. Each block is made short enough so that the signal may be 
considered stationary within the segment. The Fourier transform is then applied to each 
segment given by; 

S(Tj) = , (3.6) 

where w*(r-T) denotes the sliding window, * represents the conjugation and 
S(T,f) displays the evolution of the signal’s frequency over time. The STFT represents a 
compromise between the time-and frequency-based views of a signal. It provides 
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information about time and spectral behavior of a signal. However one can only obtain 
this information with a fixed precision, where the precision is determined by the window. 
Many different window functions, w(t), may be used. The choice will effect the 
characteristics of STFT. Once a window function has been chosen its shape and length 
will determine the resolution in time and in frequency. As a result of the uncertainty 
principle, the time resolution (Af), and the frequency resolution (A/) of a given signal 
are inversely related. Their product has lower bound of l/4;r, which is achieved by the 
Gaussian window [9]. This produces a trade off of time resolution for frequency 
resolution and vice versa. Since the choice of window will fix (Ar)and (A/) for the 
entire time axis, the STFT partitions the time-frequency plane into a uniform grid. The 
window can not simultaneously provide good time resolution (requiring short windows) 
and good frequency resolution (requiring long windows). The performance (i-.e., 
frequency resolution) of a window is governed by its functional form [22]. 

B. WAVELET ANALYSIS 

1. Introduction 

A wavelet is defined as an oscillatary function of time or space [10]. The term 
wavelet comes from the French word ondelette, which means small wave. It has its 
energy concentrated in time, which allows the analysis of transient, non-stationary, or 
time-varying phenomena. We will take a wavelet transform and use it in the series 
expansion of signals similar to the Fourier series approach. 

2. The Continuous Time Wavelet Transform (CTWT) 

The Fourier transform of g(t) is defined as; 

= . (3.7) 
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This is the integral over all time of the signal g(t) multiplied by a complex exponential. 
Similarly, the continuous time wavelet transform (CTWT) is defined as the integral over 
all time of the signal multiplied by scaled, shifted versions of a wavelet function '?(/): 


C(t,<i)=4=L«)W— 
Va L I a 


dt. 


(3.8) 


where 'F(r) is the wavelet function. The parameter t denotes translation in time, and the 


scale factor a denotes dilation or compression in time. The factor 1 /^fa normalizes the 
energy of the CTWT and * denotes conjugation. 

One advantage of the wavelet analysis is that it allows the selection of large 
number of basis functions, contrary to being restricted to sinusoids as in the Fourier 
analysis. Two important characteristics of the wavelets are that; the wavelet function 
'F(r) is of finite duration, the wavelet function 'F(0has a zero mean and is zero at the 
end points. The first and second characteristic requires that the basis function oscillates 
about zero, and gives rise to the name wavelet or small wave [10]. 

The time resolution and frequency resolution of the CTWT is controlled by the 
scale factor a (Equation 3.8). Low scales (small values of a) correspond to high 
frequency wavelets and provide good time resolution. High scales (large values of a) 
correspond to low frequency wavelets with poor time resolution but good frequency 
resolution. 

A second advantage of the wavelet analysis is the multi-resolution capability it 
provides in time-scale plane. The time-frequency mapping of the STFT and the CWT is 
shown in Figure 3.1 
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Figure 3.1: (a) Frequency-Time plane for STFT, (b) Scale-Time plane for CWT. 


The STFT produces a uniform grid with a constant time (Ar) and frequency resolution 
(A/), while the CTWT has a time and frequency resolution that depends on the scale. 
Note that the CTWT has a time resolution that improves at higher frequency while 
frequency resolution degrades. 

3. The Discrete Wavelet Transform (DWT) 

Although the discretized continuous wavelet transform enables the computation 
of the continuous wavelet transform by computers, it is not a true discrete transform [6]. 
As a matter of fact, the wavelet coefficients are simply a sampled version of the CWT, 
and the information it provides is highly redundant, as far as the reconstruction of the 
signal is concerned. This redundancy, on the other hand, requires a significant amount of 
computation time and resources. 
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The DWT is sufficient for most practical applications for the reconstruction of the 
signal. The DWT provides sufficient information, and can offer a significant reduction in 
the computation time. It is considerably easier to implement than the continuous wavelet 
transform and obtained by restricting the scale and time parameters of the CWT to 
discrete values. The DWT of a discrete signal g(n) is defined by 

C(a,^) = , (3.9) 

n ^ja y a ^ 

where a, b, and n are the discrete versions of a, t, and t of Equation 3.8 respectively. The 
scaling factor is further restricted to; 

a = ai, J=0,l,...log2(iV). (3.10) 

The choice of agWill govern the accuracy of the signal reconstruction via the inverse 
transform. It is popular to choose a^-l, since it permits the implementation of fast 

algorithms. Setting a = 2^ produces octave bands called dyadic scales. As the scale level 
is increased from J to J+1, the analysis wavelet is stretched in the time domain by a factor 
of two. Hence the DWT output has better frequency resolution and less precise time 
resolution as the scale number increases. 

If the time shifting parameter b is restricted to k2 ^, where k is an integer, this 
version of the DWT is known as the decimated DWT and can be written as: 

C(2\k2') = y-^g(nyV(2-'n-k)- (3.11) 

n=l ■yCl 

where 7 =0,1, 2,-", logjCA), ^ =1, 2, ••-, iV2“^ , and N is the length of the signal g(n). 

The term k2^ in the argument of DWT, indicates that C(a, b) is decimated by a factor of 
two at each successive scale J by retaining only the even points. 


a. 


Subhand Coding and Multiresolution Analysis 


The time-scale (frequency) representation of the signal is obtained by 
using digital filtering techniques. The CWT is computed by changing the scale of the 
analysis window, shifting the window in time, multiplying by the signal, and integrating 
over all times. In the discrete case, filters of different cutoff frequencies are used to 
analyze the signal at different scales. The signal is passed through a series of high pass 
filters, and passed through a series of low pass filters. Each one of the filter is followed 
by a two-to-one decimator. 

The DWT analyzes the signal at different frequency bands with different 
resolutions by decomposing the signal into a coarse and a detail component at each scale. 
The DWT employs two sets of functions, a scaling function and a wavelet function. 
These can be associated with a lowpass and a highpass filter, respectively. The 
decomposition of the signal into different frequency bands is obtained by successive 
highpass and lowpass filtering of the signal followed by the decimation operation. 

Figure 3.2 illustrates this procedure. Here g[n] is the original signal to be 
decomposed, and d[n] and h[n] are the lowpass and highpass filters respectively. The 
bandwidth of the signal at every level is marked on the figure as BW. The original signal 
is first passed through a halfband highpass filter h[n] and a lowpass filter d[n]. After 
filtering, half of the samples can be eliminated according to the Nyquist’s rule. This is 
because the output has a high frequency oi nil radians instead of n. The signal can be 
subsampled, by simply discarding every other sample. This decomposition halves the 
time resolution, since the output is characterized by half the number of samples compared 
to the original signal. However this operation doubles the frequency resolution, since the 


26 



frequency band of the signal now spans only half the input frequency band. The above 
procedure, which is known as the subband coding, is repeated. 

C. THE EFFECTINESS OF WAVELET ANALYSIS 

Wavelet expansions and wavelet transforms have been shown to be very efficient 
and effective in analyzing a very wide class of signals and phenomena [10]. 

1. The wavelet expansion allows a more accurate time description and 
identification of signal characteristics. A Fourier coefficient represents a component that 
lasts for the integration time of the transform and, therefore, temporary events must be 
described by a phase characteristic that allows cancellation or reinforcement over large 
time periods. A wavelet expansion coefficient represents a component that is itself local 
and is easier to interpret. The wavelet expansion may allow a separation of components 
of a signal that overlap in time or frequency. 

2. Wavelets are adjustable and adaptable. Because there is not just one 
wavelet, they can be designed to fit individual applications. They are ideal for adaptive 
systems that adjust themselves to accommodate the signal. 

3. The generation of wavelets and the calculation of the discrete wavelet 
transform are well suited for digital implementation. 
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Figure 3.2: The Subband coding algorithm (All bandwidth refers to the original sampling 
rate). , 
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IV. TIME DIFFERENCE OF ARRIVAL (TDOA) ESTIMATION 


In this chapter, we will discuss how to estimate the TDOA between two signals. 
The TDOA can be employed to find the position of a GSM emitter. Figure 4.1 shows a 
typical configuration; one emitter and a pair of receivers. The signals are shown in their 
discrete time form. 



The correlation function can be used to estimate the TDOA. In the next section, 
we will discuss about how to calculate the correlation function and to determine the 
TDOA from it. 
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A. CORRELATION FUNCTION 


Frequently we would like to know the association between two signals, that is, 
how one signal is related to the other. Correlation of signals is often encountered in radar 
and sonar processing, digital communications, and other areas of science and 
engineering. 

To be specific, let us suppose that we have situation, as shown in Figure 4.1. The 
signal at receiver one is denoted by x(k), while y(k} represents the time shifted version of 
x{k) at receiver 2. With additive noise, x(k) and y{k) can be modeled as: 
x{k)=s{k)+n^{k) (4.1) 

y{k)=as{k-D)+n^{k) it = 0,1, •••, N-1 , (4.2) 

where s(k) is the unknown source signal, nj(fc)and n 2 {k)zxt additive noises at the 
receivers, D is the difference in arrival times at the receivers, a is an attenuation 
coefficient, and N is the number of samples in each snap shot received at the two 
receivers. 

The most widely accepted method for obtaining TDOA (£> in Equation 4.2) uses 

the cross correlation method. Expectation of x(k) and y(k) leads to 

R^(j) = E{x(k)y(k-T)] 

= £:{[s(ifc) +n^ik)][as(k-D-T)+n2ik- t)]} 

= E\as{k) s{k-D-r) + s(k)n2 ik-T) + as(k-D- r)n^ (k) +n^ (k)n2 (A: - t)} . 

Since the noise and signal are independent, and the noise has zero mean, all the cross 
terms in the expectation equal to zero. Also the noises are independent of each other, so 
the expectation of the n^{k) and n 2 (it -T)is also zero. The theoretical cross correlation of 
x(k) and y(k) becomes 
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R^{T)=aR^{D + t) . (4.3) 

Figure 4.2 shows the circuit for cross correlation estimation and the block 


diagram for the discrete time cross correlator. The cross correlation approach requires 
that receivers share a precise time reference. The performance of the TDOA estimates 
can be improved by increasing the summation interval. Once the cross correlation 
function is computed, the value of r which maximizes Equation 4.3 is used as the 
estimate of the TDOA. 



Figure 4.2: (a) Circuit diagram for measuring the cross correlation function of x(k) and 
y(k) (b) Block diagram for the discrete time cross correlator. 
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Figure 4.3 displays the equivalent fast correlation method. A cross spectral density 
estimate is obtained in the frequency domain, and the cross correlation estimate is 
obtained via an inverse Fourier transform. The fast correlation method requires 
appropriate zero padding prior to taking the Fourier transform of x{k} and y(k). 

B. SIMULATION 

Simulations were carried out to evaluate the performance of the cross correlation 
method in performing the TDOA estimate of signals airiving at two separate receivers. 

The /'* error is given by 
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Figure 4.3; Fast cross correlation method using fast Fourier transforms. 


e, =Z)-D, i = l,2,-,N , (4.4) 

where N is the number of realizations. The mean square error (MSE) is given by 

= > < 4 . 5 ) 

Ntr' 

where e,. is the error for the i'* realization. In this part of the experiment the 13-bit Barker 
code ([+ + + + + -- + + - + - +]) was used, because it has good correlation properties. A 
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200 noise realizations were used and the delay D for this part was fixed at a value of 25 
time samples. We simulated two modulation scenarios. 

i. The symbol + was modulated by a sinusoid of one period over 32 time 
samples (the sampling frequency fs=l/32) while the symbol - phase 
shifted the sinusoid by 180 degree relative to the + symbol. 

ii. The symbol + was modulated by a sinusoid of one period over 8 time 
samples (the sampling frequency fs=l/8) while the symbol - phase shifted 
the sinusoid by 180 degree relative to the + symbol. 

Figure 4.4 plots the mean square error versus to the signal-to-noise ratio (SNR) 
for case (i) and case (ii). The mean square error depends on the sampling rate. We note 
that, as the number of samples increases relative to the bit period, the mean square error 
decreases. 
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(a) 



Figure 4.4: (a) The mean square error in estimating the TDOA using cross correlation of 
time domain signals for case i (b) for case ii. 
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V. WAVELET DENOISING 


Wavelet decompositions have many applications in signal processing. One 
important one is noise reduction. Each transform coefficient represents a measure of the 
correlation between the signal and the basis function. Large coefficients represent good 
correlation, conversely small coefficients represent poor correlation. The main idea in 
denoising is to retain the coefficients that preserve the signal while removing those that 
represent noise. 

Two main properties of the wavelet transform assist in separating the noise 
coefficients from the rest. The first property is that, by choosing the basis function to 
match the signal, the resulting decompositions will contain relatively few coefficients. 
The second property is that, for a Gaussian noise input, the transform coefficients will 
remain Gaussian. The orthogonal wavelet transform is a linear operation, which will 
transform Gaussian noise into Gaussian noise [18]. 

The wavelet shrinkage algorithm introduced by Donoho and Johnstone [18, 19, 
20] keeps only the significant coefficients, representing the signal based on non-linear 
thresholding. It discards the coefficients that fall below a given magnitude. Wavelet 
denoising consists of three steps. 

1. Taking the wavelet transform of the input signal. 

2. Selecting a threshold and suppressing the noisy coefficients by applying a 
non-linear thresholding technique. 

3. Performing the inverse wavelet transform based on the modified 
coefficients to reconstruct the signal. 
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A. 


CALCULATING A THRESHOLD VALUE 


The term threshold refers to a constant that is computed as the cut-off value to 
separate the coefficients that will be retained from those that will be suppressed. The 
noisy data can be expressed, 

x(k)=s(k)+crn(k) it =0,1, 2,-", V , (5.1) 

where s(k) is the input signal to be recovered from the noisy data, n(k) is zero mean, unit 
variance additive white Gaussian noise (AWGN), N is the length of the signal, tris the 
standard deviation of the AWGN. Since the wavelet transform is a linear operator, the 
wavelet coefficients will be in the form of Equation 5.1 

Algorithms computing the threshold value T require estimation of <7. Five 
methods of computing thresholds are described below. 

1. Stein’s Unbiased Risk Estimator (SURE) Threshold ) 

This method of threshold calculation was proposed by Donoho and Johnstone 
[20]. The thresholding is adaptive: A threshold level is assigned to each dyadic resolution 
level by the principle of minimizing the Stein Unbiased Estimate of Risk (SURE) [21] for 
threshold estimates. The SURE threshold is smoothness adaptive. If the unknown signal 
contains jumps, the reconstruction does also. If the unknown signal has a smooth 
segment, the reconstruction is as smooth as the mother wavelet will allow. This statistical 
procedure calculates the estimated mean square error (risk) for a range of threshold 
values, and selects T^y with the resulting minimum risk. 

2. Sqtwolog Threshold () 

Sqtwolog threshold uses a fixed form threshold yielding minimax performance 
multiplied by a small factor proportional to log(length(signal)) [7] 
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=^j2log(length(signal)) 


(5.2) 


3. Heursure Threshold () 

Heursure threshold is mixture of SURE and sqtwolog threshold methods [7]. 

4. Minimaxi Threshold () 

Minimaxi threshold uses a fixed threshold chosen to yield mimimax performance 
for mean square error against an ideal procedure [7]. The minimax principle is used in 
statistics in order to design estimators. It is designed to select estimators that minimize 
the worst case error occurring in the set. 

5. Wo-So-Ching Threshold () 

This threshold method was proposed by Wo, So, and Ching. It will be discussed 
in Chapter VI. 

The threshold techniques, 1 through 4, are part of Donoho’s method. Details are 
discussed in Chapters VI and Vn. Donoho’s method used at the low SNR’s did not result 
in improvement. Threshold technique 5 is used successfully and details are discussed in 
Chapter VI and Vn. 

B. THRESHOLDING METHODS (SHRINKAGE) 

Once a threshold value is established, a number of methods exist to apply the 
threshold to suppress or modify the coefficients of decomposition. We examined three 
thresholding methods. 

1. Hyperbolic Thresholding 

Hyperbolic thresholding was proposed by Wong and will be discussed in Chapter 
VI. 
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2. Hard Thresholding 

Hard thresholding can be described as the usual process of setting to zero the 
elements whose absolute values are lower than the threshold [7]. For notational 
convenience, we define x(k) and n{k) as vectors. 


X=[40),x(l), x(2),-,x(iV-l)]" 
iV=[n(0), n(l),n(2),-,n(Ar-l)]" 


Let W be NxN a wavelet transform matrix. In a vector form, the transformed output C is 
related to the input vector X by C=W X , where 

C = [c(7,0, ; = i = , (5.4) 

and J = log 2 (iV). The indices j and i represent the scale and the position in each scale, 
while c(-l,0) denotes the remaining low-pass filtered coefficients. 

The non-linear hard threshold can be defined 



c{i,j) 

0 


, for |c(i,;)|>r 
; otherwise . 


(5.5) 


As seen in Figure 5.1, hard thresholding of the transformed coefficients retains the 
coefficients that exceed the threshold value, while all others are set to zero. 

3. Soft Thresholding 

As we see in Figure 5.2, soft thresholding is an extension of hard thresholding. It 
first sets to zero the elements whose absolute value is lower than the threshold, and then 
shrinks the remaining non-zero coefficients by the threshold amount [7]. 

The non-linear soft thresholding can be defined 

[sign (c(i, 7))lc(i, ;)| - f] ; for |c(i, ;)| ^ T 


c{i,j) = 


0 


otherwise . 


(5.6) 
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The advantage of this method is that the results are not as sensitive to the precise value of 
the threshold T, as in the “keep or kill” strategy of the hard thresholding. The 
disadvantage of this method is that the general shape of the signal might be slightly 
affected since the even the large coefficients are modified when using this scheme. 

The thresholding techniques 2 and 3 are part of Donoho’s method, which did not 
lead to improvement. Thresholding method 1, hyperbolic thresholding, was used 
exclusive in our work. Details can be found in Chapter VI and VII. 


Original signal Hard thresholded signal 



Figure 5.1: Hard thresholding applied to a data segment. 
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Figure 5.2: Soft thresholding applied to a data segment. 






VI. IMPROVED TDOA ESTIMATION 


In this chapter, we will discuss the performance of several denoising schemes 
designed to improve TDOA estimation. Seven methods are presented and evaluated. 
These tend to reduce the effective noise at the receivers. The methods are: wavelet 
denoising based on Donoho’s method [7], wavelet denoising using the Wo-So-Ching 
threshold. [15], wavelet denoising using hyperbolic shrinkage [16], wavelet denoising 
using median filtering [14], a modified approximate maximum-likelihood delay 
estimation based in part on [17], denoising based on the fourth order moment, and a time 
varying technique. Figure 6.1 is a generic figure for all methods addressed in this chapter. 

A. WAVELET DENOISING BASED ON DONOHO’S METHOD 

The wavelet denoising as proposed by Donoho [7, 18, 19, and 20], was discussed 
in Chapter V. This method fails at low SNR’s and the results can be found in Chapter 
VII. 

B. WAVELET DENOISING USING THE WO-SO-CHING THRESHOLD 

Wavelet denoising is applied to the noisy signals, which are obtained at the two 
spatially separated sensors. Prior to cross correlation, each one of the sensor outputs is 
denoised according to the Wo-So-Ching thresholding rule to increase the effective input 
signal-to-noise ratio (SNR). The proposed system consists of two sub-units (Wavelet 
denoising and Correlator) as depicted in Figure 6.1. 

Wavelet denoising (WD) is applied to each received signal to recover the 
corresponding source waveform. The restored signals are cross correlated. The TDOA 
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estimate is given by the argument at which the cross correlation function attains its 
maximum value. 



Figure 6.1: System block diagram for TDOA estimate using wavelet denoising. 


The wavelet denoising technique consists of three steps; taking the wavelet 
transform of the received signals, thresholding to modify the wavelet coefficients, and 
performing the inverse wavelet transform of the modified coefficients. For notational 
convenience, we define the x(k),y(k),s(k),n^{k),andn2ik) sequences in vector form as 

X = [x(0), x(l), x(2), ■■■,x(N-l)V 
F=[y(0), y(l), yi2),-,yiN-l)V 

5 = k0), r(l), si 2 ), -, s(N-l)V (6.1) 

A^i = [«i (0), n, (1), n, (2), •••,«! (A^ -1)] ^ 

N2=[n2 (0), (1), (2), • • •, (iV -1)] ^ 

The start time “0” denotes the first data point in a given block of data. Let W be a NxN 
orthonormal wavelet transform matrix. The transformed output C is related to the input 
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vector X by C=W X . Since W is a linear transform matrix, C can be decomposed into 
C=C^ , where 

C,=[c^(j,i), j = i = 0,l,2,--,2"-'] (6.2) 

7 = -1,0,1,--,7;i = 0,1,2,---,2''’] , (6.3) 

are the wavelet transform coefficients of the source signal vector S and the noise vector 
Ni^ (k=l, 2) respectively. The indices i and j represent the scale and the position in each 

scale, respectively, while Cj(-1,0) and (-1,0) denote the low-pass filtered coefficients. 

Note that, C is still a Gaussian vector. The main idea of the signal restoration 
using wavelet denoising is to adapt each to make its value close to Cg(j,i) so that 

a good approximation of s(k) can be obtained after taking the inverse wavelet transform. 

The Wo-So-Ching threshold is derived according to the Neyman-Pearson 
criterion as it is used in hypothesis testing [15]. This criterion is stated as follows: 

Let c be a Gaussian random variable with known variance . Let a test be conducted 
with the following hypotheses 

H,:E{c} = iu, (6.4) 

and, 

H,:E{c}^Mo (6-5) 

We will denote the acceptance of hypothesis //q and as Dq and Dj, respectively. The 
type n error PiD^ / i/j) will be minimized for a given P(Dj / /fg) when the threshold A is 
selected as follows 

\c-JilQ\<A = ^^2<J.elf'\l-a) , ( 6 . 6 ) 
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where 


c is the Gaussian random variable with variance , 
(Iq is the expectation of c under hypothesis, 

is the Wo-So-Ching threshold, 
a is the type I error that is «= f’(D, and 


erf{v)=^\e‘^dt 

0 


(6.7) 


This denoising method discards the individual elements in Q that are too small in 
magnitude. The wavelet coefficient c(j,i) is regarded as totally due to noise if 
^c{ j,i)\< A, = 'j2(y erf ~\l-a) . (6.8) 

As a result, c, (j,i) will be modified and is given by 


c,(aO= 


fc(7,0 , if \cU,i)\^A, 


0 


otherwise 


(6.9) 


The denoising process is applied to both x(k) and y(k). The restored signals of s(k) 
mda s(k - D) are denoted by 5(fc) and a sik - D ), respectively. The TDOA estimate is 
given by the location of peak of the cross correlation function of modified x(k) and y(k} 
coefficients. 

C. WAVELET DENOISING USING HYPERBOLIC SHRINKAGE 

As mentioned earlier, denoising of signals corrupted by white noise uses three 
steps; the wavelet transform, shrinkage, and the inverse wavelet transform. Shrinkage 
reduces the value of each coefficient in magnitude by an amount related to the threshold 
value. The amount of the shrinkage is determined by the shrinkage method. In this 
section the hyperbolic method is used. 
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Hyperbolic shrinkage [16] is defined as: 




sign[cii, j)]^|cii,jy 
0 


|c(z, 7 )|<i ; 


( 6 . 10 ) 


where c(i,j) is the detail function of received signal, y^-is the Wo-So-Ching threshold ( as 
given in Equation 6.8), and c(z, j) is the modified detail function after shrinkage. The 
coefficients are set to the square root of the difference of the squares of the values as long 
as the absolute value of the coefficient is greater than the threshold. If the absolute value 
of the coefficient is less than the threshold then the coefficient is set to zero. Hyperbolic 
shrinkage can give good performance for a wide variety of signals. Other shrinkage 
techniques may offer better performance for specific input signals and noise conditions. 
The proposed denosing system is shown in Figure 6.1 
D. WAVELET DENOISING USING THE MEDIAN FILTER 

In this method, a median filter (filter length 3) is applied to the first detail function 
of the wavelet transforms. The signals are reconstructed using the modified wavelet 
coefficients. The denoising technique is depicted in Figure 6.1. 

The median filter is a one dimensional filtering technique. It is a non-linear 
technique that applies a sliding window to a sequence [14]. The median filter replaces the 
center point of the window with the median value of all the points contained in the 
window. The length of the window is very important. For example, for a stationary signal 
such as a sinusoid, a longer window length is better. But if the signal is a non-stationary, 
a short window length is better. This is an important drawback when using the median 
filter if one does not have a priori information about the source signal. 
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E. MODIFIED APPROXIMATE MAXIMUM-LIKELIHOOD DELAY 
ESTIMATION 

An approximate maximum likelihood (AML) algorithm is proposed by Y.T. 
Chan, H.C. So, and P.C. Ching to estimate the TDOA of signals [17]. The general idea in 
this method is to attenuate the frequency bands where the noise is strong and to enhance 
the frequency bands where the signal is strong. To do this, we can use the pre-filters as 
shown in Figure 6.2. The pre-filters H^{f) and // 2 (/)are chosen as follows [17] 


GssU) 


h,(/)H2(/)=— 


1+ 


GJf) , GJf) 


( 6 . 11 ) 




Where G„(/), G„_„^(/),and denote the auto-power spectra of s(k), riiik)and 


iij(it)respectively. This choice of known as maximum likelihood (ML) 

weighting. 

When we multiply the denominator and nominator in Equation 6.11 by 
Gn,n, (f)G„^„, (/) and integrate the power spectral densities, we obtain the formula for the 
weight for each subband. This tends to enhance frequency bands where the signal is 
strong. 

The weights are given by 


w. = 


sdi 


^n2di ^sdi ^ 


/ j^2 


^2 




( 6 . 12 ) 


and are used as indicated in Figure 6.3. 
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Figure 6.2: A generalized cross correlator. 


Since the noise density is flat and the filter gain of the wavelet transform is two, 
the noise power will essentially remain same after passing through each high and low 
pass filter. The noise power at each subband can be estimated using the formulas below 



(6.13) 


(6.14) 

=argmax^^(T) 

(6.15) 

A 7 


(6.16) 

^2 

^n2di ” ^/i2 

(6.17) 
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The variance are the noise powers at each receiver, while is the signal power 

at receiver, represent the noise powers at the/"' subband after the wavelet 

transform. The signal power at each subband “/” can be estimated by 


a-2 _ 


/ 1 N-l 


1 /V-l 




/t=0 


_ ^2 I = 1 2 • • • 7 • 


(6.18) 


where N is the length of detail function d,., and “/” denotes the scale. 


The denoising procedure, shown in Figure 6.3, consists of the three steps. These 
steps are discrete wavelet decomposition, scaling of each subband sequence, and an 
inverse wavelet transform. The modified subband sequences are obtained by weighting 
each subband by 

= d,. ; / = l,2,-,7 . (6.19) 

The modified AML method is applied to both channels. The modified subband 
sequences and a^are then combined to form the modified AML denoised 

signals. The MATLAB codes for modified approximate maximum-likelihood technique 
are presented in the Appendix. 

F. WAVELET DENOISING BASED ON THE FOURTH ORDER MOMENT 


1. White Noise 

White noise is the term applied to any zero mean random process whose power 
spectral density is constant. Since the power spectral density function is constant, the 
correlation function for white noise is an impulse. In other words, the samples of white 
noise process are uncorrelated. The correlation function of white noise process M[n] has 
the form 
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1 > «o ] = {n]S[n^ - «o ] 


(6.20) 


Usually a white noise process is taken to a mean a stationary white noise process so that 

K\.n^(TlS[l] (6.21) 

and 




( 6 . 22 ) 


If Mis a random vector of N consecutive samples from white Gaussian process, 
the probability density function is given by 




1 


20’o 


N-\ 

=n 


1 


JifJ: 


(6.23) 


«=o 

The samples of a Gaussian white noise process are uncorrelated and hence independent. 



Figure 6.3; Block diagram for scaling of each subband. 
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2. 


Moments of White Noise Process 


The definitions of moments for a stationary random process are; 


=£{»[«]}= (6.24) 

[/, ]= E\u{n\u[n +]}= [/j ] (6.25) 

Af f ’ [/,, Zj ]= E\u[n]u[n + /, ]«[« + /j ]} (6.26) 

+ + + • (6-27) 


The first two moments are equal to the mean and correlation function defined earlier. The 
fourth moment is 

^ [Zj, Z2» ^3 ] = + ^2 + Z3 ]} 

= E{u.[n]u[n + Z, ]}£■{«[« +12 ]«[« + Z3 ]} ^g) 

+ E\u[n]u[n + Zj ]}£{«[« + Zj ]M[n + Zj ]} 

+ £{«[«]«[« + Z 3 ]}£{«[« + Z, ]«[« + Zj]} . 


[0,0,0] =£’{M[n]M[n]}E{M[n]M[n]} 

+ £'{M[n]M[n]}£'{M[n]M[n]} 
+ E{M[n]M[n]}£{M[n]M[22]} 


(6.29) 


3. The Wavelet Transform of White Noise 

Since the wavelet transform is a linear operation, the Gaussian process remains 
Gaussian after passing through low and high pass filters and hence preserves the high 
order moments properties. 

4. Fourth Order Moment of the Received Signal 


The received signals are of the form 
x(k) = s(k)+ni{k) 
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(6.30) 




y{k) = a s(k - £>)+Wj (^) 


(6.31) 


where n^(k)and n2(k)aie, statistically independent, zero mean Gaussian random 

variables. The fourth order moment estimation of the received signal is, 

[0,0,0] =£:{x:" [it]} 

=E{[s[k] + n,[k]Y} 

=£{5n^]}+4£:t'A:K[^]}+6£{5'[/:K[^]} 

+ 4£:{s[it]nf [it]}+ e { i * [it]} 


5. 


Mean and Standard Deviation of the Fourth Order Moment 


If we call z = MY’ , the mean of the z becomes, 


JV ,-=o 1=0 


(6.33) 


The second moment of the z is. 


f ^ f 1 W-l N-\ 

[iV /=o 7=0 

1 


iV “If iV 




N N N 


(6.34) 


The variance of the z is, 

a\ 

='^(W+96)-9<7.» 

= 56,,. , 

N ‘ 


(6.35) 
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If N is large the, variance is close to zero. We accounted for the variance in choosing a 
threshold, see Equation 6.42. Experimentally, this amounted to a change from 3cr^^ to 

3.1a". 

'*1 

6. Denoising Based on the Fourth Order Moment 

Since the wavelet transform is linear, each detail function in Figure 6.3 is in the 
same form as Equation 6.30 and 6.31. The fourth order moment of detail function which 
contains signal should be greater than 3 <t^ , where denotes the noise power at 

subband c?,.. Using this property we can eliminate the wavelet coefficients which 

represent noise and hold the ones which represent the signal. 

We can estimate the noise power using the formulas below 
<7^ = argmax R„(t) (6.36) 

r 

<=K(0)-ct^ (6.37) 

. (6.38) 

The noise power after passing through the first high pass filter is 

a;,,=G.cr„M2 7 = 1,2 , (6.39) 


where G is the gain of the high and low pass filters. If we generalize this formula, the 
noise power of any detail function is 


And we assume that filter gain is 2 then 

7 = 1,2 


(6.40) 


^2 _ — 2 
a, j =cr. 


(6.41) 
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Using the fourth order moment property we can modify the each detail function in Figure 
6.3 as given by 


d^ = 


k iif E{l‘}iOa‘,±<Tl) 


0 ; otherwise 


(6.42) 


After modifying all detail functions, we perform the inverse wavelet transform to 
reconstruct the denoised signal. The MATLAB codes for wavelet denosing based on the 
fourth order moment technique are presented in the Appendix. 


G. TIME VARYING TECHNIQUE 

The modified AML delay estimation method is based on the idea of attenuating 
the spectral components where the signal is weak. This method works well if the subband 
information is noise only or signal only for the duration of the segment. 

Since some signals are time varying or transient, the subbands could contain 
signal only for parts of the segment. We modified this method by allowing for some 
variations over time. In particular, we divided each detail output into two time segments 
and applied the equations of section E to each segment. This approach allowed a time 
varying gain (i.e weights for the first half of the segment, w^j weights for the second 

half of the segment). The MATLAB codes for time varying modified approximate 
maximum-likelihood technique (time varying MAML) are presented in the Appendix. 
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VII. SIGNAL DESCRIPTION AND SIMULATION RESULTS 


A. SIGNAL DESCRIPTION 

We used two different types of test signals, a generic signal and a GSM signal. 
These two signals are described in this section. 

1. Generic Signals 

In this part of the simulation we used the 13 bit Barker code sequence as the 
message to create four different test signals. The generic signal representing one code bit 
is of the form sin(2 * f*nlN) where N=32, n=0,1,31. Signals A, B, C, and D use 
value of/of 1, 4, 8, and 12 respectively. The generic signals are plotted in Figure 7.1. 
The MATLAB codes for generic signal generation are presented in Appendix. 


(a) Test Signal A (b) Test Signal B 



Figure 7.1: (a) Generic signal A, (b) Generic Signal B, (c) Generic signal C, (d) Generic 
signal D. 
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2. GSM Signal 

The GSM signal properties are detailed in Chapter 11. According to these 
properties we simulated the GSM signal. The I and Q channel of the GSM base band 
signal are shown in Figure 7.2. The MATLAB codes for GSM signal generation are 
presented in Appendix [4]. 


(a) I channel 



Figure 7.2: (a) I channel of GSM signal (b) Q channel of GSM signal. 
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3. Delay DeHnition 

For the generic signal simulations the delay was randomly chosen between -30 to 
+30. For the GSM signal simulations the delay was randomly chosen between -150 to 
+150. The realizations of delays are shown in Figure 7.3 in a histogram fashion. Both 
randomizations used the same seed number. 


(a) 

401 - 1 - 1 - 1 ^- 1 -,- f 



-401-1-•-1-1-1- i -1-1_I_I 

0 20 40 60 80 100 120 140 160 180 200 


(b) 

200 \ - 1 -j- 1 - 1 -r- 



-200 ^-*-'_'_'_'_‘'_*_"_ * 

0 20 40 60 80 100 120 140 160 180 200 

realization number 


Figure 7.3: (a) Delay for generic signals, (b) Delay for the GSM signals. 


B. SIMULATION RESULTS 

We used the mean square error (MSB) as the criteria of goodness (i.e. low MSB 
denotes a small error and hence good localization). The MSB is defined in Bquation 4.5. 
To quantify the global improvement, we compute the sum of the MSB’s obtained from a 
given method and compare it to the total error for the non-denoised TDOA estimation. 
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The improvement is given in percentage. The total error is obtained by summing the 
MSE over all SNR’s of interest (i.e., -3 to +5 dB). 

1. Donoho’s Method 

Donoho proposed different types of thresholds and thresholding techniques. This 
is explained in Chapter V. All of this methods were tried on the generic signals. The 
mean square error using wavelet denoising based on Donoho’s method did not improve 
relative to the non-denoised version at the low SNR values of interest. The mean square 
error obtained from the correlation of the time domain raw signals is much smaller than 
the ones obtained using Donoho’s method. 

2. Wavelet Denoising Using the Wo-So-Ching Threshold 

a. Simulation Results for the Generic Signals 

As seen in Figure 7.4 (a), there is an slight improvement in the mean 
square error using Wo-So-Ching threshold method for generic signal A. As the carrier 
frequency is increased (Figure 7.4 (b) and Figure 7.5), the mean square error also 
increase. As long as the number of samples per binary bit is high there is an improvement 
using this denoising method, but the advantage disappears as the carrier frequency 
increases. The results are sununarized in Table 7.1. 

b. Simulation Results for the GSM Signals 

As seen in Figure 7.6, there is about 28% improvement in the total mean 
square error using Wo-So-Ching thresholding method for GSM signal. 
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(a) 



Figure 7.4: (a) MSE versus SNR with and without denoising using the Wo-So-Ching 


thresholding method for generic signal A. (b) MSE for generic signal B. 


(a) 



Figure 7.5: (a) MSE versus SNR with and without denoising using the Wo-So-Ching 


threshold method for generic signal C. (b) MSE for generic signal D. 
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SNR 

Xcorr 

(A) 

Wo-So-Ching 

threshold 

(A) 

Xcorr 

(B) 

Wo-So-Ching 

threshold 

(B) 

5 

0.385 

0.25 

0 

0 

4 

0.415 

0.275 

0 

0 

3 

0.485 

0.425 

0 

0 

2 

0.715 

0.535 

0 

0 

1 

0.795 

0.78 

0.32 

0.32 

0 

0.805 . 

0.725 

0 

0.96 

-1 

1.145 

1.02 

2.89 

3.84 

-2 

1.62 

1.34 

4.505 

6.49 

-3 

2.085 

1.815 

8.38 

12.36 


Table 7.1: (a) MSE versus SNR for with and without denoising using the Wo-So-Ching 


threshold method for generic signal A and B. 


SNR 

Xcorr 

(C) 

Wo-So-Ching 

threshold 

(C) 

Xcorr 

(D) 

Wo-So-Ching 

threshold 

(D) 

5 

0 

0 

0 

0 

4 

0 

0.16 

0 

0 

3 

0 

0.56 

0 

0 

2 

0.08 

0.64 

0 

0 

1 

1.2 

3.52 

0 

0 

0 

2.72 

5.84 

0.32 

1.92 

-1 

3.52 

8.56 

1.6 

5.12 

-2 

5.92 

16.56 

5.095 

9.7 

-3 

8.88 

28.08 

5.665 

15.62 


Table 7.1: (b) MSE versus SNR for with and without denoising using the Wo-So-Ching 


threshold method for generic signal C and D. 
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SNR 

Figure 7.6: MSE versus SNR with and without denoising using the Wo-So-Ching 
threshold method for the GSM signal. 

3. Wavelet Denoising Using the Hyperbolic Shrinkage 
a. Simulation Results for the Generic Signals 
As seen in Figure 7.7, there is an improvement in the mean square error 
using hyperbolic shrinkage method for generic signal A and B. Again as the carrier 
frequency is increased (Figure 7.8), the mean square error also increase. The results are 
summarized in Table 7.2. If we compare the Table 7.1 and Table 7.2, we can conclude 
that the hyperbolic shrinkage method gives better result than the Wo-So-Ching threshold 
method. 


61 






b. Simulation Results for GSM Signals 

As seen in Figure 7.9, there is about 48% improvement in the total mean 
square error using Hyperbolic Shrinkage method for GSM signal. If we compare the 
Figure 7.6 and 7.9 we can conclude that hyperbolic shrinkage method gives better result 
for GSM signal than Wo-So-Ching thresholding method. 


(a) 



SNR 


Figure 7.7: (a) MSE versus SNR with and without denoising using the Hyperbolic 
Shrinkage method for generic signal A. (b) MSE for generic signal B. 
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(b) 



SNR 


Figure 7.8: (a) MSB versus SNR with and without denoising using the Hyperbolic 


Shrinkage method for generic signal C. (b) MSB for generic signal D. 


SNR 

Xcorr 

(A) 

Hyperbolic 

shrinkage 

(A) 

Xcorr 

(B) 

Hyperbolic 

shrinkage 

(B) 

5 

0.375 

0.2000 

0 

0 

4 

0.53 

0.3200 

0 

0.08 

3 

0.52 

0.4100 

0 

0.08 

2 

0.625 

0.5400 

0 

0.16 

1 

0.96 

0.7350 

0.32 

0.24 

0 

1.07 

0.8200 

0 

1.12 

-1 

1.345 

1.0050 

2.89 

2.8 

-2 

1.36 

1.2100 

4.505 

3.6 

-3 

1.895 

1.3300 

8.38 

7.84 

-4 

2.425 

2.2400 

12.025 

10.465 

-5 

2.57 

3.5350 

38.24 

17.185 

-6 

87.41 

14.5850 

268.01 

190.41 


Table 7.2: (a) MSB versus SNR for with and without denoising using Hyperbolic 
shrinkage method for generic signals A and B. 
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SNR 

Xcorr 

(C) 

Hyperbolic 

shrinkage 

(C) 

Xcorr 

(D) 

Hyperbolic 

shrinkage 

(D) 

5 

0 

0.1200 

0 

0 

4 

0 

0.3200 

0 

0.24 

3 

0.08 

0.4800 

0 

0.32 

2 

0.16 

1.0000 

0 

1.76 

1 

1.04 

2.0400 

0 

2.8 

0 

1.6 

3.2200 

0 

4.015 

-1 

4.08 

4.3800 

1.645 

7.625 

-2 

4.08 

5.0800 

0.96 

12.735 

-3 

9.76 

10.6450 

8.76 

87.35 

-4 

14.8 

16.6750 

15.385 

29.055 

-5 

24.96 

20.1500 

22.97 

680.8 

-6 

118.8 

163.7250 

57.635 

1278.2 


Table 12\ (b) MSE versus SNR for with and without denoising using Hyperbolic 
shrinkage method for generic signals C and D. 



Figure 7.9: MSE versus SNR with and without denoising using the Hyperbolic Shrinkage 
method for the GSM signal. 
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4. Wavelet Denoising Using the Median Filter 

a. Simulation Results for the Generic Signals 

As seen in Figure 7.10, there is an improvement in the mean square error 
using median filtering on generic signals A and B. Again as the carrier frequency is 
increased (Figure 7.11), the mean square error also increase. The results are summarized 
in Table 7.3. If we compare the Table 7.1, 7.2 and 7.3, we can conclude that median 
filtering has the best result for the generic signal B but for the other signals the 
Hyperbolic shrinkage method provides better results. 


(a) 



Figure 7.10: (a) MSE versus SNR with and without denoising using the median filtering 
method for generic signal A. (b) MSE for generic signal B. 
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Figure 7.11: (a) MSE versus SNR with and without denoising using the median filtering 
method for generic signal C. (b) MSE for generic signal D. 


SNR 

Xcorr 

(A) 

Median filtering 
(A) 

Xcorr 

(B) 

Median filtering 
(B) 

5 

0.375 

0.26 

0 

0 

4 

0.53 

0.325 

0 

0 

3 

0.52 

0.49 

0 

0 

2 

0.625 

0.635 

0 

0 

1 

0.96 

0.645 

0.32 

0 

0 

1.07 

0.675 

0 

0.005 

-1 

1.345 

0.865 

2.89 

0.005 

-2 

1.36 

1.235 

4.505 

2.24 

-3 

1.895 

1.81 

8.38 

2.895 

-4 

2.425 

2.125 

12.025 

8.855 

-5 

2.57 

2.855 

38.24 

19.965 

-6 

87.41 

28.35 

268.01 

200.06 


Table 7.3: (a) MSE versus SNR for with and without denoising using the median filtering 
method for generic signals A and B. 
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SNR 

Xcorr 

(C) 

Median filtering 
(C) 

Xcorr 

(D) 

Median filtering 
(D) 

5 

0 

10.52 

0 

67.555 

4 

0 

11.32 

0 

264.57 

3 

0.08 

10.84 

0 

301.49 

2 

0.16 

11.32 

0 

963.45 

1 

1.04 

13.88 

0 

973.63 

0 

1.6 

14.11 

0 

3251.7 

-1 

4.08 

21.825 

1.645 

3593.1 

-2 

4.08 

23.425 

0.96 

5698.7 

-3 

9.76 

34.415 

8.76 

5088.8 

-4 

14.8 

178.33 

15.385 

9083.3 

-5 

24.96 

385.47 

22.97 

9397.2 

-6 

118.8 

1217.5 

57.635 

9916 


Table 7.3: (b) MSE versus SNR for with and without denoising using the median filtering 
method for generic signals C and D. 

b. Simulation Results for GSM Signals 

As shown in Figure 7.12, there is about 42% improvement in the total 
mean square error using the median filtering method for the GSM signal. If we compare 
the Figure 7.6, 7.9 and 7.12 we can conclude that hyperbolic shrinkage method gives the 
best result for the GSM signal. 

5. Modified Approximate Maximum-Likelihood Delay Estimation 
a. Simulation Results for the Generic Test Signals 
As seen in Figure 7.13 and 7.14, there is a significant improvement in the 
mean square error for all generic signals when using the modified AML estimation 
method. The results are summarized in Table 7.4. If we Compare the Table 7.1, 7.2, 7.3 
and 7.4, we can conclude that the modified AML estimation method gives the best result. 
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Figure 7.12: MSE versus SNR with and without denoising using the median filtering 
method for the GSM signal. 


(a) 



SNR 

Figure 7.13: (a) MSE versus SNR with and without denoising using the modified AML 
estimation method for generic signal A. (b) MSE for generic signal B. 
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(b) 



Figure 7.14: (a) MSE versus SNR with and without denoising using the modified AML 
estimation method for generic signal C. (b) MSE for generic signal D. 


SNR 

Xcorr 

(A) 

Modified AML 
(A) 

Xcorr 

(B) 

Modified AML 
(B) 

5 

0.375 

0.07 

0 

0 

4 

0.53 

0.085 

0 

0 

3 

0.52 

0.17 

0 

0 

2 

0.625 

0.185 

0 

0 

1 

0.96 

0.275 

0.32 

0 

0 

1.07 

0.305 

0 

0 

-1 

1.345 

0.42 

2.89 

0.325 

-2 

1.36 

0.54 

4.505 

2.255 

-3 • 

1.895 

0.695 

8.38 

4.43 

-4 

2.425 

0.815 

12.025 

7.995 

-5 

2.57 

1.26 

38.24 

13.085 

-6 

87.41 

1.635 

268.01 

61.795 


Table 7.4: (a) MSE versus SNR for with and without denoising the modified AML 
estimation method for generic signals A and B. 
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SNR 

Xcorr 

(C) 

Modified AML 
(C) 

Xcorr 

(D) 

Modified AML 
(D) 


0 

0 

0 

0 


0 

0 

0 

0 


0.08 

0 

0 

0 


0.16 

0.08 

0 

0 

1 

1.04 

0.4 

0 

0 

0 

1.6 

1.44 

0 

0 

-1 

4.08 

2.96 

1.645 

0.125 


4.08 

4.685 

0.96 

1.28 


9.76 

7.525 

8.76 

2.925 


14.8 

10.45 

15.385 

10.75 


24.96 

20.395 

22.97 

13.35 


118.8 

2670.8 

57.635 

208.15 


Table 7.4: (b) MSE versus SNR for with and without denoising using the modified AML 


estimation method for generic signals C and D. 


b. Simulation Results for GSM Signals 

As seen in Figure 7.15, there is about 79% improvement in the total mean 
square error using the modified AML estimation method relative to the undenoised 
version. Comparing Figure 7.6, 7.9, 7.12 and 7.15 we can conclude that the modified 
AML estimation method gives the best result for the GSM signal. 

6. Wavelet Denoising Based on the Fourth Order Moment 
a. Simulation Results for the Generic Signals 
As seen in Figure 7.16 and 7.17, there is a significant improvement in the 
mean square error using wavelet denoising based on the fourth order moment method for 
all generic signals. The results are summarized in Table 7.5. If we compare this method 
with the other methods, we can conclude that wavelet denoising based on the fourth order 
moment is the second best method. 
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SNR 


Figure 7.15: MSE versus SNR with and without denoising using the modified AML 


estimation method for the GSM signal. 



(b) 



Figure 7.16: (a) MSE versus SNR with and without denoising using the wavelet 


denoising based on the fourth order moment method method for generic signal A. (b) 


MSE for generic signal B. 
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(a) 



Figure 7.17; (a) MSB versus SNR with and without denoising using the wavelet 
denoising based on the fourth order moment method method for generic signal C. (b) 
MSB for generic signal D. 


SNR 

Xcorr 

(A) 

Fourth order 
(A) 

Xcorr 

(B) 

Fourth order 
(B) 

5 

0.375 

0.08 

0 

0 

4 

0.53 

0.11 

0 

0 

3 

0.52 

0.16 

0 

0 

2 

0.625 

0.215 

0 

0 

1 

0.96 

0.35 

0.32 

0 

0 

1.07 

0.47 

0 

0.32 

-1 

1.345 

0.58 

2.89 

0.64 

-2 

1.36 

0.82 

4.505 

1.6 

-3 

1.895 

0.92 

8.38 

4.73 

-4 

2.425 

1.455 

12.025 

5.775 

-5 

2.57 

1.625 

38.24 

10.58 

-6 

87.41 

2.54 

268.01 

21.2 


Table 7.5: (a) MSB versus SNR for with and without denoising wavelet denoising based 
on the fourth order moment method for generic signals A and B. 
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SNR 

Xcorr 

(C) 

Fourth order 
(C) 

Xcorr 

(D) 

Fourth order 
(D) 

5 

0 

0 

0 

0 

4 

0 

0 

0 

0 

3 

0.08 

0 

0 

0 

2 

0.16 

0.24 

0 

0 

1 

1.04 

0.72 

0 

0 

0 

1.6 

1.12 

0 

0 

-1 

4.08 

3.84 

1.645 

0.125 

-2 

4.08 

5.52 

0.96 

1.92 

-3 

9.76 

8.49 

8.76 

4.845 

-4 

14.8 

12.375 

15.385 

12.19 

-5 

24.96 

20.86 

22.97 

16.595 

-6 

118.8 

26.99 

57.635 

24.85 


Table 7.5; (b) MSE’versus SNR for with and without denoising using wavelet demising 
based on the fourth order moment method for generic signals C and D. 
b. Simulation Results for GSM Signals 

As seen in Figure 7.18, there is about 63% improvement in the mean 
square error using wavelet denoising based on the fourth order moment method. If we 
compare Figure 7.6, 7.9, 7.12, 7.15 and 7.18 we can conclude that wavelet denoising 
based on the fourth order moment ranks second in performance for the GSM signal. 

7. Time Varying Technique 

As we explained in Chapter VI, we can address signals that are non-stationary 
(i.e., the signal does not stay in a given frequency band for the length of the segment or 
has modulation characteristics). In principle, in any given subband we try to find several 
weights. This allows us to keep the information for the frequency band when the signal is 
strong. We tried this method on the GSM signal and compared it with the modified AML 
estimation method. Figure 7.19 allows a comparison of the modified AML method with 
the time varying version in which the segments are partitioned into two parts for every 
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scale. The improvement in the total mean square error using the time varying modified 
AML method is about 81%. There is no drastic improvement relative to the modified 
AML, but we believe that if the signal has a time varying property, the new method based 
on the time varying technique will give better results. 



Figure 7.18: MSE versus SNR with and without denoising using the wavelet demising 
based on the fourth order moment method for the GSM signal. 
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mean square error 



SNR 


Figure 7.19: MSE versus SNR with and without denoising using the time varying 
modified AML and modified AML technique on the GSM signal. 
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VIII. CONCLUSION AND FUTURE WORK 


A. CONCLUSION 

The accuracy of algorithms used, to localize wireless communication units, was 
studied. Chapter I reviewed why there is a need to locate cellular transmitters and the 
existing localization systems. The time difference of arrival (TDOA) method was used to 
locate emitters. Estimation of the TDOA was based on the cross correlation function. To 
increase the accuracy of the TDOA estimation, the noise in the received data was 
reduced. The wavelet transform was employed to minimize the noise, and several 
denoising methods were examined. 

The MSB’s for all methods for the GSM signal were plotted in Figure 8.1. The 
modified AML delay estimation method provided the best result, wavelet denoising using 
the fourth order moment method ranked second. All methods provided better results than 
the correlation of the raw time domain signals. The time varying technique outperformed 
the modified AML method for SNR values greater or equal to -2 dB. We believe that if 
the signal has time varying properties, the time varying technique will be the superior 
method. 

The probability of no-error for the GSM signal was plotted in Figure 8.2 and 8.3. 
The Wo-So-Ching, the hyperbolic shrinkage and the median filtering techniques, see 
Figure 8.2, gave worse results than the correlation of the raw time domain signals. In 
these techniques, a threshold value was calculated and then according to the threshold, 
each coefficient was modified. These three techniques can modify the coefficients, which 
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represent the signal in the subband. Hence we reduced the mean square error by using 
these three techniques but we also decreased the probability of no-error. 


The modified AML, the fourth order moment and the time varying AML 
techniques, see Figure 8.3, gave better results than the correlation of the raw time domain 
signals. In these techniques, we evaluated the whole subband. The subband or portion of 
it, which represented the signal was kept. The subband, which represented noise only was 
eliminated. Since there was a small chance to modify signal coefficients representing the 
signal, these three techniques improved the probability of no-error. The time varying 
AML method gave the highest probability of no-error. 



SNR 

Figure 8.1: Plot of all denoising methods for the GSM signal. 
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Figure 8.2: Plot of the probability of no-error for the Wo-So-Ching threshold, Hyperbolic 


Shrinkage and Median filtering techniques. 


B. FUTURE WORK 

1. We worked exclusively with the base band GSM signal. A follow on 
study should address GSM signals at the RF or IF level, to evaluate if there is potential 
improvement of the MSE. 

2. In our experiments, only white Gaussian noise was simulated. Fading • 
signals should be included in a follow on study. 
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SNR 


Figure 8.3: Plot of the probability of no-error for the Time varying AML, Fourth order 
moment and Median filtering techniques. 
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APPENDIX 


MATLAB CODES 


A. SIGNAL GENERATION 
1. Generic Signal 

%*5lfSr*x*Vr***lV*****x*ifc***Vf***Sr*x***x*Vf**ilrVc**Vf****Vr**Vr^>r*x***x*Vf*x:*rVr*x*iir*:r* 

%****:Sr***TSr***************************************:Jr***-************^***** 

% u_inod: This function creates generic signals 
% f=l for Test signal A 

% f=4 for Test signal B 

% f=8 for Test signal C 

% f=12 for Test signal D 

% 

% SYNTAX: y=u_inod{data) 

% 

% INPUT: data = Digital data stream 
% 

% OUTPUT: y = PSK modulated data 
% 

% SUB^FUNC: None 
% Written by Unal Aktas 

%************:Jr*********************^lr*****iJr******************^Tfr**-******* 

%************T!f****Vc*****l»fT<f***-)ir***^***^***lif***TSr**T»:*****:<c********X*l^***** 

function y=u_mod(data); 

n=0:31; % NUMBER OF SAMPLES 

one=sin{2*pi*n*f*/32); 

% f=l for Test signal A 
% f=4 for Test signal B 
% f=8 for Test signal C 
% f=12 for Test signal D 

zero=-l.*one; 
mod__data= [ ] ; 
for i=l:length(data) 
if data(i)==l 

mod_data=[mod_data one]; 
else 

mod_data=[mod_data zero]; 
end 

end 

plot(mod_data);title('modulated signal') 
y=mod_data; 
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2. GSM Signal 


^Vf*******lir***********************************'fr***^**Vf****iir***<r***’!ir^**^if* 

% GSM_SET: This script initializes the values. 

% 

% SYNTAX; gsni_set 

% 

% INPUT: 

% OUTPUT: 

% 

% 

% 

% 

% 


None 

Configuration variables created in memory, these are; 
Tb(= 3.692e-6) 

BT(= 0.3) 

OSR{= 4) 

SEED(= 931316785) 

INIT_L(= 260) 


% SUB,„FUN: None 

% Written by Jan H. Milc)celsen / Arne Norre Ekstrom 


Tb = 3.692e-6; 

BT = 0.3; 

OSR = 4; 

% INITIALIZE THE RANDOM NUMBER GENERATOR. 

% BY USING THE SAME SEED VALUE IN EVERY SIMULATION, WE GET THE SAME 
% SIMULATION DATA, AND THUS SIMULATION RESULTS MAY BE REPRODUCED. 

% 

SEED = 931316785; 
rand('seed',SEED); 

% THE NUMBER OF BITS GENERATED BY THE DATA GENERATOR. (data_gen) 

% 

INIT_L = 114; 

% SETUP THE TRAINING SEQUENCE USED FOR BUILDING BURSTS 
% 

TRAINING = tOOlOOlOlllOOOOlOOOlOOlOlll]; 

% CONSTRUCT THE MSK MAPPED TRAINING SEQUENCE USING TPvAINING. 

% 

T_SEQ = T_SEQ_gen(TRAINING); 
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^x*^*********************7r******:****'*********'Ar***x***********x*-**w***-r* 

% burst_( 


% 

% 

% 

% 

% 


This function generates a bit sequence representing 
a general GSM information burst. Included are tail 
and Ctrl bits, data bits and a training sequence. 

The GSM burst contains a total of 148 bits accounted 
for in the following burststructure (GSM 05.05) 


% 

[ TAIL 

DATA 

CTRL 

TRAINING 

CTRL 

DATA 

TAIL 

] 

% 

[ 3 

57 

1 

1 26 1 

1 

57 

3 

] 


% 

% SYNTAX: 
% 

% INPUT: 


[TAIL] = [000] 

[CTRL] = [0] or [1] here [1] 

[TRAINING] is passed to the function 

tx_burst = burst_g(tx_data, TRAINING) 

tx_data: The burst data. 

TRAINING: The Training sequence which is to be used. 


% OUTPUT: tx_burst: A complete 148 bits long GSM normal burst binary 

% sequence 


% GLOBAL: 

% 

% SUB_FUNC: None 

% Written by Jan H. Mikkelsen / A.rne Norre Ekstrom 


function tx_burst = burst_g(tx_data, TRAINING) 

TAIL =[00 0]; 

CTRL = [1]; 

% COMBINE THE BURST BIT SEQUENCE 
% 

tx_burst = [TAIL tx_data(1:57) CTRL TRAINING CTRL tx^data(58:114) 
TAIL] ; 


83 


ir i( ic ir rk 'k ic :k if'-'k i( ie ir •k ic i( 'k if ic -k ir 'k it -k ic if -k ie ir ■■k it ir ic ir ir ir i( i( ic Tk if if -k ir "k -k -k ie if "k -k i( •)( -ic -k ir ir it -k ic "k 'k -k i( -k -k -k 




% 

% data„_.gen: This function generates a block of random data bits. 

% 


% SYNTAX: [tx„data] = data^gen (INIT__L) 

% 

% INPUT: INIT_L: The length of the generated data vector. 

% 

% OUTPUT: tx_data: An element vector contaning the random data 

% sequence of length INIT__L. INIT_L is a variable 

% set by gsm_set. 


% 

% SUB^FUNC: None 

% Written by Jan H. Mikkelsen / Arne Norre Ekstrom 

<^kkifkkki<-kkkifkkkirkkkifkkki(kkkifkkkifkkki(kkkitkkkitkkkitkkki(kkkiekkk'kkkkickkkitkkk 


function 
tx_data = 


[tx_data] = data_gen{INIT_L); 
round(rand(1,INIT_L)); 
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% GSM_MOD: 


% 

% 

% 

% 

% 


This MatLab code generates a GSM normal burst by- 
combining tail, Ctrl, and training sequence bits with 
two bloks of random data bits. 

The data bits are convolutional encoded according 
to the GSM recommendations 

The burst sequence is differential encoded and then 
subsequently GMSK modulated to provide oversampled 
I and Q baseband representations. 


% 

% SYNTAX: [ tx_burst , I , Q ] = gsm_mod (Tb,OSR,BT,tx_data,TRAINING) 


% 

% INPUT: Tb: Bit time, set by gsm_set.m 

% OSR: Oversampling ratio (fs/rb), set by gsm_set.m 

% BT: Bandwidth Bittime product, set by gsm_set.m 

% tx_data: The contents of the datafields in the burst to be 

% transmitted. Datafield one comes first. 

% TRAINING: The Training sequence which is to be inserted in the 

% burst. 

% 


% OUTPUT: 

% tx_burst: The entire transmitted burst before differential 
% precoding. 

% I: Inphase part of modulated burst. 

% Q: Quadrature part of modulated burst. 

% Written by Jan H. Mikkelsen / Arne Norre Ekstrom 

%********************************************************************** 

%****»*******************«***»*******«***************************«***»* 


function [ tx_burst , I , Q ] = gsm_mod{Tb,OSR,BT,tx_data,TRAINING) 

tx_burst = burst_g(tx_data,TRAINING); 

% DIFFERENTIAL ENCODING. 

% 

burst = diff_enc(tx_burst); 

% GMSK MODULATION OF BURST SEQUENCE 
% 

[I,Q] = gmsk_mod(burst,Tb,OSR,BT); 
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% diff_enc: 


% 

% 

% 

% 


This function accepts a GSM burst bit sequence and 
performs a differential encoding of the sequence. The 
encoding is according to the GSM 05.05 recommendations 

Input: D(i) 

Output: A(i) 


D'(I) = D(i) (+) D(i-l) ; (+) = XOR 

D{i), D'(i) = {0,1} 

A(i) = 1 - 2*D'(i) 

A{i) = {-1,1} 

SYNTAX: [diff_enc_data] = diff_enc(burst) 

INPUT: burst A binary, (0,1), bit sequence 

OUTPUT: diff_enc_data A differential encoded, (+1,-1), version 

of the input burst sequence 


% 


SUB_FUNC: None 

Written by Jan H. Mikkelsen / Arne Norre Ekstrom 


function DIFF_ENC_DATA = diff_enc(BURST) 


L = length(BURST); 

% INTERMEDIATE VECTORS FOR DATA PROCESSING 


d_hat = zeros(1,L); 
alpha = zeros(l,L); 

% DIFFERENTIAL ENCODING ACCORDING TO GSM 05.05 
% AN INFINITE SEQUENCE OF I'ENS ARE ASSUMED TO 
% PRECEED THE ACTUAL BURST 
% 

data = [1 BURST]; 
for n = 1+1:(L+1), 

d_hat(n-l) = xor( data(n),data(n-1) ); 
end 

alpha = 1 - 2.*d_hat; 

% PREPARING DATA FOR OUTPUT 
% 

DIFF_ENC_DATA = alpha; 
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%*i^r:Sf*x*Vf**:frVr***<**x*n;r*******T!r*******Vf*******^** *********-A *x*****Vr******* 
^************************************************;i:***********^**Vc****** 
% 

% gmsk^mod: This function accepts a GSM burst bit sequence and 

% performs a GMSK modulation of the sequence. The 

% modulation is according to the GSM 05.05 recommendations 

% 


% SYNTAX: [i.q] = gmsk_mod(burst,Tb,osr,BT) 


% INPUT: 
% 

% 

% 

% 

0.3) 

% 


burst A differential encoded bit sequence (-1,+!) 

Tb Bit duration (GSM: Tb = 3.6926-6 Sec.) 

osr Simulation oversample ratio, osr determines the 

number of simulation steps per information bit 
BT The bandwidth/bit duration product (GSM: BT = 


% OUTPUT: 
% 

% 


i,q In-phase (i) and quadrature-phase (q) baseband 

representation of the GMSK modulated input burst 
sequence 


% SUB_FUNC: ph_g.m This sub-function is required in generating the 

% frequency and phase pulse functions. 

% Written by Jan H, Mikkelsen / Arne Norre Ekstrom 

^****************************************************************i^r***** 


^********************************************************************** 


function [I,Q] = gmsk_mod(BURST,Tb,OSR,BT) 

[g/Q] = ph_g(Tb,OSR,BT); 

% PREPARE VECTOR FOR DATA PROCESSING 
% 

bits = length(BURST); 

f_res = zeros(1,(bits+2)*OSR); 

% GENERATE RESULTING FREQUENCY PULSE SEQUENCE 
% 

for n = l:bits, 

f_res((n-l)*OSR+l:(n+2)*OSR) = f_res((n-1)*OSR+l:(n+2)*OSR) + 
BURST(n).*g; 
end 

% CALCULATE RESULTING PHASE FUNCTION 
% 

theta = pi*cumsum(f_res) ; 

% PREPARE DATA FOR OUTPUT 
% 

I = cos(theta); 

Q = sin(theta); 
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%******»*************************************************************** 

********************************************************************* 

% 

% PH_G: This function calculates the frequency and phase functions 

% required for the OMSK modulation. The functions are 

% generated according to the GSM 05.05 recommendations 


% SYNTAX: [g_fun, g_fun] = ph_g (Tb, osr,BT) 


% INPUT: 

% 

% 

% 

% 

% OUTPUT: 
% 

% 


Tb Bit duration (GSM: Tb = 3.692e-6 Sec.) 

osr Simulation oversample ratio, osr determines the 

number of simulation steps per information bit 
BT The bandwidth/bit duration product {GSM; BT = 0.3) 

g_fun, q_fun Vectors contaning frequency and phase 

function outputs when evaluated at osr*tb 


% SUB^FUNC: None 

% VJritten by Jan H. Mikkelsen / Arne Norre Ekstrom 

^★Tlr****************-*--*-**!^****-*************'**** ********* ^’'f***'***'^***'^***"** 


function [G_FUN, Q_FUN] = ph_g(Tb,OSR,BT) 
Ts = Tb/OSR; 

% PREPARING VECTORS FOR DATA PROCESSING 


PTV = -2*Tb:Ts:2*Tb; 
RTV = -Tb/2:Ts:Tb/2-Ts; 


% GENER.ATE GAUSSIAN SHAPED PULSE 
% 

sigma = sqrt(log(2))/(2*pi*BT); 

gauss = (l/(sqrt(2*pi)*sigma*Tb))*exp{-PTV.'-2/(2*sigma^2*Tb''2)); 

% GENERATE RECTANGULAR PULSE 
% 

rect = 1/(2*Tb)*ones(size(RTV)); 

% CALCULATE RESULTING FREQUENCY PULSE 
% 

G_TEMP = conv(gauss,rect); 

% TRUNCATING THE FUNCTION TO 3xTb 
% 

G = G_TEMP(OSR+l:4*OSR) ; 

% TRUNCATION IMPLIES THAT INTEGRATING THE FREQUENCY PULSE 
% FUNCTION WILL NOT EQUAL 0.5, HENCE THE RE-NORMALIZATION 
% 

G_FUN = (G-G(l))./(2*sum(G-G{l))); 

% CALCULATE RESULTING PHASE PULSE 
% 

Q_FUN = cumsum(G_FUN) ; 
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B. DENOISING TECHNIQUES 


1. Modified AML 


^****:Sf***************TSr***********^***********-r*******Vr***Tir****^/r******-!fr* 

% amll: Approximate maximum-likelihood delay estimation 
% via orthogonal wavelet transform. In this function, 

% we assumed that noise is Gaussian noise and it has a 

% flat freq response. We modify each detail function 

% by multipliying modified AML coefficient based on 

% the signal and noise powers 

% 

% SYNTAX: y=amll(xn,yn,delay) 

% 

% INPUT: xn = Received signal from first receiver 
% yn = Received signal from second receiver 

% delay = True TDOA betv^een two received signals 
% 

% OUTPUT: 

% y = Error between true TDOA and estimated TDOA 

% 

% SUB^FLINC: None 
% Written by Unal Aktas 


function y=amll(xn,yn,delay); 

xyn=xcorr (xn,yn, 'biased') ; 

[sigmas b]=max(xyn); 
rx=xcorr(xn,'biased'); 
maxx=rx(length(xn)); 
ry=xcorr(yn,'biased'); 
maxy=3ry( length (yn) ) ; 
sigmanl=maxx-sigmas ; 
sigman2 =maxy-sigmas; 

nx=floor (log2 (length (xn) ) ) ,* 
ny=floor(log2(length(yn))); 

[cx lx]=wavedec(xn,nx,'db4'); 

[cy ly] =wavedec (yn,ny,' db4') ; 

dxc=[]; 
for i=l:nx 

d=detcoef(cx,lx,i); 
dl=length(d); 
sigmad= (1/dl) *sum(d. '^2) ; 
sigmasd=sigmad-sigmanl; 
if sigmasd<=0 
wd=0; 
else 

wd=sigmasd/ (sigmanl'*^sigman2+sigmasd* (sigmanl+sigman2) ) ; 
end 
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dc=wd*d; 
dxc=[dc dxc]; 

end 

a=appcoef(cx,lx,'db4',nx); 
al=length(a); 
sigmaa= (1/al) *s\im(a. ^2) ; 
siginasa=sigmaa-siginanl; 
if siginasa<=0 
wa=0; 
else 

wa=sigmasa/ (sigmanl*sigman2+siginasa* (sigmanl+sigman2) ) ; 
end 

ac=wa*a; 
dxc=[ac dxc]; 

xd=waverec(dxc,lx,'db4'); 

dyc=[]; 
for i=l:ny 

dy=detcoef(cy,ly,i); 
dyl=length(dy); 
sigmady= (1/dyl) *suin(dy. ^"2) ; 
siginasdy=sigmady-siginanl ; 
if sigmasdy<=0; 

wdy=0; 

else 

wdy=sigmasdy/ (sigmanl*sigman2+siginasdy* (sigmanl+sigman2) ) 
end 

dcy=wdy*dy; 
dyc= [dcy dye] ; 

end 

ay=appcoef(cy,ly,'db4',ny); 
ayl=length(ay); 
siginaay= (1/ayl) *sum(ay. '"2) ; 
siginasay=siginaay“Siginanl ; 
if siginasay<=0 
way=0; 
else 

way=sigmasay/ (sigmanl*sigman2+sigmasay* (sigmanl+sigman2 ) ) ; 

end 

acy=way*ay; 
dyc=[acy dye]; 

yd=waverec(dye,ly/'db4'); 

rxyd=xcorr(xd,yd); 

[a5 b5]=max(rxyd); 

er5=delay“(b5-1024) ; 
y=er5; 
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% 

% tezS: This is a test program for modified AML technique. In this 
% program, we used GSM signal. 

% 

% SYNTAX: tez5 
% 

% INPUT: None 
% 

% OUTPUT: Mean sqaure error versus SNR 
% 

% SUB__FUN: amll.m 
% Written by Unal Aktas 

C^ieiciciv-kieici^ieiticificic-k'k-kieic-k'k-k'k'k'k-k'k'k-k'kir-k'kic'k'k^ic-k'k-irie'k^r-k-k-k'k-kit-k'k-ificirif'kiriT-k-k’kir'k'kic'h'k'kif 

clear all 
gsm_set; 

data=data_gen{INIT_L); % this creates a binary data 
[tx_burst/1,Q]=gsm_mod{Tb,OSR,BT,data,TRAINING) ; 

s=I+j*Q; 
sl=length(s); 

pow= (1/sl) *sum(abs (s) . ^2) ; 

K=200 % NUMBER OF REALIZATIONS 

rand('seed',40); 
f=150*rand{l,K); 
delay=floor(f); 

delay (l:K/2) =-delay(l:K/2) ; % DELAY IS BET&VEEN -150 TO +150 

n=[5 4 3 2 1 0 -1 -2 ~3 -4 -5 -6]; 

SNR=10.^(n./lO) ; 

for k=l:length(SNR) 

oran=SNR(k) 
for i=l:K 

x=[zeros(1,200) s zeros(1,224)]; 

y=[zeros(1,200+delay(i)) s zeros(1,224-delay(i))]; 
randn{'state',2*(i+j)); 

noil_real=sqrt(pow/(2*oran))*randn(l,1024); 
randn('state’,3 *(i+j)); 

noil_imag=sqrt(pow/(2*oran))*randn(1,1024); 
randn('state‘,4*(i+j)); 

noi2_real=sqrt(pow/(2*oran))*randn(l,1024) ; 
randn('state',5*(i+j)); 

noi2_iinag=sqrt (pow/ (2*oran) ) *randn (1,1024) ; 

noil=noil_real +j *noil_imag; 
noi2=noi2_real+j *noi2__imag; 
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xn=x+noil; % x ,+ noise 
yn=y+noi2; % y + noise 

% TDOA CALCULATION WITHOUT DENOISING 

xy=xcorr(xn,yn); 

[al bl]=max(real(xy)); 

erl(i)=delay(i)-(bl-1024); 

%modified AML 

el=amll(real(xn),real(yn),delay(i)); 
e2=amll (iinag(xn) ,imag(yn) , delay (i) ) ; 

er8(i)=(el+e2)/2; 


end 

error8a(k) = (1/length (er8)) *sum(er8 .''2) ; 
errorla(k)=(l/length(erl))*sum(erl .^2); 

H8a(k,:)=er8; 

Hla(k,;)=erl; 


end 

figure(6) 

k=[5 4 3 2 1 0 -1 -2 -3 -4]; 

plot(k,errorla(1:10),'o',k,errorSa(1:10) , 'x',k,errorla(1:10),k,.. 
error8a(l:10)) 

legend('xcorr without denoising','modified AML') 

save errorSa; 
save errorla; 

save Hla; 
save H8a; 


92 



2 , 


Time Varying MAML 


^**^*Vr**i^r:«r***:rr***:*:*******x***:Jr***********************^***:A:* ********.A- 
^rfr**i«r****T!r*******iir*-*-*«*******Ttr***********:rr***<r***^***T^********* ******** 

% ami2 : We modified the MAML technique by dividing each 
% detail function into two segments. Different 

% coefficients for each segment are computed. 

% 

% SYNTAX: y=aml2(xn,yn,delay) 

% 

% INPUT: xn = Received signal from first receiver 

% yn = Received signal from second receiver 

% delay = True TDOA between xn and yn 

% OUTPUT: y = Error between true TDOA and estimated TDOA 

% 

% SUB^FUNC: None 
% Written by Unal Aktas 

%********************************************************** **;^:***y.***** 
%********************************************************************** 


function y=aml2(xn,yn,delay); 

xyn=xcorr{xn,yn,'biased'); 

[sigmas b]=max(xyn); 
rx=xcorr(xn,'biased'); 
maxx=rx(length(xn)); 
ry=xcorr(yn,'biased') ; 
maxy=ry(length(yn)); 
sigmanl=maxx-sigmas; 
s igman2 =maxy-s i gmas ; 

nx=floor(log2(length(xn))); 
ny=floor(log2(length(yn))); 

[cx lx]=wavedec(xn,nx,'db4'); 

[cy ly]=wavedec(yn,ny,'db4'); 

dxc= [ ] ; 
for i=l;nx 

d=detcoef(cx,lx,i); 
dl=length(d); 

Nsl=128/2'"(i-l) ; % THE LENGTH OF THE SUBBLOCK 

if Nsl<=l 
Ns=dl; 
else 

Ns=Nsl; 

end 

D=ceil(dl/Ns); 
if dl<Ns*D 

dm=[d zeros (l,D*NS“dl) ] ; 

end 

for k=l:D 

p=(k-1)*Ns+l:k*Ns; 
sigmad= (1/Ns) *sum(dm(p) .'"2) ; 
sigmasd=sigmad-sigmanl; 
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if siginasd<=0 
wd=0; 
else 

wd=sigmasd/ (sigmanl*siginan2+sigmasd* (sigmanl+sigman2) ) ; 
end 

dc (p) =wd*din(p) ; 

end 

dxc= [dc (1 :dl) dxc] ; 

end 

a=appcoef(cx,lx,'db4',nx); 
al=length(a); 
siginaa= (1/al) * sum (a. ^2) ; 
siginasa=sigmaa-siginanl; 
if sigmasa<=0 
wa=0; 
else 

wa=siginasa/(sigTnanl*siginan2+siginasa* (siginanl+sigman2 ) ) ; 

end 

ac=wa*a; 
dxc=[ac dxc]; 

xd=waverec (dxc, lx, ' db4 ') ; 

dyc=[]; 
for i=l;ny 

dy=detcoef (cy, ly, i) ; 
dyl=length(dy); 

Nsl=128/2"(i-l); 
if Nsl<=l 
Ns=dyl; 
else 

Ns=Nsl; 

end 

D=ceil(dyl/Ns); 
if dyl<Ns*D 

diny=[dy zeros (l,D*Ns) ] ; 

end 

for k=l:D 

p=(k-l)*Ns+l:k*Ns; 
sigmady= (l/Ns) *suin(diny(p) .'" 2 ) ; 
sigmasdy=siginady-sigmanl ; 
if sigmasdy<=0; 

wdy=0; 

else 

wdy=siginasdy/ (siginanl*sigman2+siginasdy* {siginanl+siginan2) ) 
end 

dcy (p) =wdy*diny (p) ; 

end 

dyc= [dcy(l:dyl) dye] ; 

end 

ay=appcoef (cy, ly, 'db4' ,ny) ; 

ayl=length(ay); 

sigmaay=(1/ayl)*sum(ay.^2); 
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sigmasay=siginaay-sigmanl ; 
if sigmasay<=0 
way=0; 
else 

way=sigmasay/ {sigmanl*sigman2+sigmasay* (siginanl+sigman2) ) 

end 

acy=way*ay; 
dyc=[acy dye]; 

yd=waverec(dye,ly,'db4'); 

rxyd=xeorr(xd,yd); 

[a5 b5 ] =max {rxyd) ; 

er5=delay*-(b5-1024) ; 
y=er5; 
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% tez5t : This is a test program for time varying MAML technique. 

% In this program, GSM signal is used 

% 

% SYNTAX: tezBt 
% 

% INPUT; None 
% 

% OUTPUT: Mean square error versus SNR 
% 

% SUB_FUNC: aml2.m 
% Written by Unal Aktas 

^ififififififificifififififificificifirifitifififififififififififififificififkifififificififkififificififififififififififififififitificifif 
if if if if if if if if if if if it if if if if if if if if if if if if if if if if if if if if if if if if if if if if if if if if if if if if it if if if if if if k it if if k if if if if it if if k if 

clear all 
gsm_set; 

data=data_gen(INIT_L); 

[ tx_burst, I, Q] =gsm_mod (Tb, OSR, BT, data, TRAINING) ; 

s=I+j*Q; 
sl=length(s); 

pow=(1/sl)*sum(abs(s).^2); 

K=200 % NUMBER OF REALIZATIONS 

rand('seed',40); 
f = 150*rand(l,K) ; 
delay=floor(f); 

delay(l:K/2)=~delay (l:K/2) ; % DELAY IS BET^'JEEN “150 TO +150 


n=t5 4 3 2 1 0 -1 -2 >3 “4 “5 -6]; 

SNR=10.'"(n./10) ; 

for k=l:length(SNR) 

oran=SNR(k) 
for i=l:K 

[zeros(1,200) s 2eros(l,224)]; 
y=[zeros(1,200+delay(i)) s zeros(1,224-delay(i))]; 

randn(‘state',2*(i+j)); 

noil_real=sqrt(pow/(2*oran))*randn(l,1024) ; 
randn('state‘,3*(i+j)); 

noil_imag=sqrt(pow/(2*oran))*randn(l,1024); 
randn('state',4*(i+j)); 

noi2_real=sqrt(pow/(2*oran))*randn(l,1024) ; 
randn('state',5*(i+j)) ; 

noi2_imag=sqrt(pow/(2*oran))*randn(1,1024) ; 

noil=noil_real+j *noil_imag; 
noi2=noi2_real+j *noi2_imag; 
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xn=x+noil; % x + noise 
yn=y+noi2; % y + noise 

% TDOA CALCULATION WITHOUT DENOISING 

xy=xcorr{xn,yn); 

[al bl]=max(real(xy)); 

erl(i)=delay(i)~(bl-1024); 

%Time varying MAML 

el=aml2 (real (xn) , real (yn) ,delay (i) ) ; 
e2=aml2 (imag{xn) , imag(yn) ,delay (i) ) ; 

er8(i)=(el+e2)/2; 


end 

errorSt (k) = (1/length(er8) ) *sum(er8 .'^2); 
errorla(k)=(l/length(erl))*sum{erl.^2); 

H8t(k,:)=er8; 

Hla(k,:)=erl; 


end 

load error8a 
figure(6) 

k=: [5 4 3 2 1 0 -1 -2 -3 -4] ; 

plot(k,errorla(1:10),'o',k,error8t(1:10),'x',k/error8a(1:10), 'd' 
k,errorla(l:10),k,error8t(1:10),k,error8a(1:10)) 

legend('xcorr without denoising','time varying ami','ami') 

save error8t; 
save errorla; 

save Hla; 
save H8t; 
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3 . 


Wavelet Denoising Based On The Fourth Order Moment 


^.****************-#r***:V***:^***Sr***iV***************Sr********-Ar*,*r^*******T^* 

% sta: Wavelet Denosing Based on The Fourth Order Moment 
% 

% SYNTAX: y=sta(xn,yn,delay) 

% 

% INPUT: xn = Received signal from first receiver 
% yu = Received signal from second receiver 

% delay = True TDO.A between two xn and yn 

% 

% OUTPUT: y = Error between true TDOA. and estimated TDOA 
% 

% SUB_FUNC: None 
% Written by Unal Aktas 

%**«***************************»*******************************«******* 

^***’************i»r************:ir*******************Tlr******T«r**-*********i«f** 

function y=sta(xn,yn,delay); 

xyn=xcorr(xn,yn,'biased'); 

[sigmas b]=max(xyn); 
rx=xcorr(xn,'biased'); 
maxx=rx{length(xn)); 
ry=xcorr(yn,'biased'); 
maxy=ry(length(yn)); 
sigmanl=maxx“sigmas; 
s igman2 =maxy ~ s igmas ; 

lamdax=3 .l*sigmanl'"2; 
lamday=3 .l*sigman2'"2; 

nx=floor(log2(length(xn))); 
ny=floor(log2(length(yn))); 

[cx lx]=wavedec(xn,nx,'db4'); 

[cy ly]=wavedec(yn,ny,'db4'); 

dxc= [ ] ; 
for i=l:nx 

d=detcoef(cx,lx,i); 
dl=length(d); 

A=(1/dl)*sum(d.^4) ; 

if A<lamdax 

dc= 2 eros(l,dl) ; 
else 

dc=d; 

end 

dxc=[dc dxc); 

end 

a=appcoef(cx,lx,'db4',nx); 
al=length(a); 
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A= (1/al) *suin(a. ^4) ; 


if A<lamdax 

ac=zeros(1,al); 
else 
ac=a; 

end 

dxc=[ac dxc]; 

xd=waverec(dxc,lx,'db4'); 

dyc= [ ] ; 
for i=l:ny 

dy=detcoef(cy,ly,i); 
dyl=length(dy); 

A=(l/dyl)*sum(dy.^4); 

if A<lamday 

dcy=zeros{1,dyl); 
else 

dcy=dy; 

end 

dyc= [dcy dye] ; 

end 

ay=appcoef (cy, ly, 'db4' ,ny) 
ayl=length(ay); 

A= (1/ayl) *siiin(ay. ^4) ; 

if A<lamday 

acy=zeros(1,ayl); 
else 

acy=ay; 

end 

dyc=[acy dye]; 

yd=waverec(dye,ly,'db4'); 

rxyd=xeorr(xd,yd); 

[a5 b5]=max(rxyd); 

er5=delay-(b5-1024); 
y=er5; 


^x**************:************************'*:***********************-* ****** 

% tezS: This is a test program for wavelet denosing 
% based on the fourth order moment technique. 

% GSM signal is used. 

% 

% SYNTAX: tezG 
% 

% INPUT: None 
% 

% OUTPUT: Mean square error versus SNR 
% 

% SUB^FUNC: sta 
% Written by Unal Aktas 

^********************************************************************** 
%****************************★***************************************** 
clear all 

gsm_set 

data=data_gen (INIT__L) ; 

[ tx_burst, I, Q] =gsm_mod(Tb, OSR,BT, data,TRAINING) ; 

s=I+j*Q; 
sl=length(s); 

pow= (1/sl)*sum(abs(s) .^2); 

K=200 

f=150*rand(l,K) ; 
delay=floor(f); 
delay(l;K/2)=-delay(l:K/2); 
n=[5 43210-1-2-3-4-5 -6]; 

SNR=10.'"{n./10) ; 

for k=l:length(SNR) 

oran=SNR(k) 
for i=l:K 

x=[zeros(1,200) s 2eros(l,224)]; 

y=[zeros(1,200+delay(i)) s zeros(1,224-delay(i))]; 
randn{'state',2*(i+j)); 

noil_real=sqrt(pow/{2*oran))*randn(1,1024) ; 
randn('state',3*(i+j)); 

noil_imag=sqrt(pow/(2*oran))*randn(1,1024) ; 
randn('state',4*(i+j)); 

noi2_real=sqrt(pow/(2*oran))*randn(1,1024) ; 
randn('state' , 5*(i + j)) ; 

noi2_imag=sqrt(pow/(2*oran))*randn{l,1024) ; 

noil=noil_real+j *noil_imag; 
noi2=noi2_real+j *noi2_imag; 
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xn=x+noil; 
yn=y+noi2; 


% TDOA CALCULATION WITHOUT DENOISING 

xy=xcorr{xn,yn); 

[al bl]=max(real(xy)); 

erl(i)=delay(i)-(bl-1024) ; 

%WAVELET DENOISING BASED ON THE FOURTH' ORDER MOMENT 

el=sta{real(xn),real(yn),delay(i)); 
e2=sta(imag(xn),imag(yn),delay (i)); 

erl0(i)=(el+e2)/2; 


. end 

errorlOa(k) = (1/length(erlO) ) *sum(erlO . ^^2) ; 
errorla(k) = {l/length{erl) ) *suin(erl. ^"2) ; 

H10a{k,:)=erlO; 

Hla{k,:)=erl; 

end 

figure(6) 

k= [5 4 3 2 1 0 -1 -2 -3 -4] ; 

plot.(k,errorla(l:10) , 'o' ,k, errorlOa (1:10) , 'x' , k, errorla (1:10) ,k,errorl0 
a{l:10)) 

legend('xcorr without denoising','Fourth Order Moment') 

save errorlOa; 
save errorla; 

save Hla; 
save HlOa; 
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